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During the last two decades, use of statistical models and methods for the analysis of dynamic phenomena 
has become widespread in various fields of natural sciences and technology. Statistical dynamics, concerned 
with the study of various random phenomena in dynamic systems, enriches the classical basic theory of 
oscillations. It also extends the possibilities for its application to the description and analysis of real response 
processes in dynamic systems, which in many cases may be of quite irregular or random nature. In fact, 
many phenomena can be described adequately only within the frame work of statistical models i.e. 
deterministic models prove to be inadequate, or at least extremely complex. Problems of this kind arise, 
frequently, in rotor dynamics, especially those involving random vibrations of bearings (rolling element and 
fluid film both), excited by random loads. Moreover, rotor bearing problems are inherently strongly nonlinear 
(Dimentberg, 1988; Bendat, 1990; Soong and Grigoriu, 1993; Lin and Cai, 1995; Ragulskis et al, 1974; 
Harris, 1984; Stolarski, 1990; Zhou and Hashimoto, 1995; Choyetal, 1992; Childs, 1993). 

Rolling element bearing are known to possess highly non-linear elastic characteristics and the complexities 
and approximate nature of the analytical determination of these properties are responsible for some of the 
unreliability in the prediction of the response and stability of a rotor system. Procedures [Ragulskis, et. al., 
1974; Harris, 1984] are available for estimation of bearing stiffness under static loading conditions in 
isolation of shaft. The existing parameter estimation techniques for overall rotor-bearing systems (e.g. 
Muszynska and Bently, 1990) involve, a controlled input excitation , that needs to be given to the bearings. 

The present study attempts for the development of a parameter estimation procedure, for the non-linear 
elastic parameters of bearings, that can be used on-line, and is based on the analysis of easily accessible 
vibration signals, picked up from the bearing caps. 



The forces setting the rotor-bearing system comprise of both - harmonic forces, due to rotor unbalance and 
random forces due to bearing surface imperfections, caused by random deviations from their standard 
theoretical design and progressive surface and subsurface deterioration. In addition excitation can be 
contributed by from random inaccuracies in the rotor-bearing-housing assembly etc. Stochastical dynamics 
methods provide an opening, to attempt the inverse problem of parameter estimation, in such a rotor-bearing 
system, governed by a non-linear equation of motion, having both deterministic and random excitation. Under 
certain engineering assumptions, the random excitation can be taken to be an ideal white noise process and 
the response of a dynamic system (linear or non-linear) can be modeled as a diffusive Markov Process or an 
approximate Markov Process. 

The structure of a Markov Process is completely determined for all future times by the distribution at some 
initial time and by a transition probability density function, which satisfies the Fokker-Planck-Kolmogorov 
(FPK) equation. The FPK equation can be solved to obtain analytical expressions for the first order 
probability distribution of the stationary response. For the response to be modeled by a Markov Process, it is 
necessary that die excitations be approximated by ideal white noises. This restriction can be, in principle, 
removed at the price of increasing the complexity of the system. Linear filters are introduced between the 
ideal white-noise excitations and the system, to produce realistic excitation spectra for the nonlinear system. 
The ideal white-noise excitations, driving the filters, guarantee that the response of the extended system is a 
Markov Process (of higher order), for which an extended FPK equation can be written. For such an 
approximate Markovian modeling stochastic averaging is applied, whereby rapid fluctuations are averaged to 
provide simpler equations for slowly fluctuating quantities. The procedure is a nontrivial extension of the 
Krylov-Bogoliubov (1937) averaging method for deterministic excitations, since it involves accounting for 
the averaged effect of a random excitation multiplied by a correlated response. 

The problem is initially attempted for balanced rotors, by considering the random forces from the bearings 
as the primary source of excitation, which are generally large enough to cause measurable level of vibrations. 
The dynamics of the rotor-bearing system is modeled as a Markov Process and Fokker-Planck equations are 
formulated. The drift and diffusion coefficients in the FPK equation are derived from the nonlinear equations 
of motion of the dynamic system (Crandall, 1958). The approach to the solution of equation is greatly 
simplified if the overall random excitation to the system, from the variety of sources, is treated as ideal white 
noise. While many engineering applications are based on this idealization, insofar as the excitation itself is 
concerned, it turns out that the response obtained through such models are quite acceptable if the time scale 



of excitation is much smaller than the time scale of the response (Lin, 1973). The time scale for the 
excitation is the correlation time, roughly defined as the length of time separation beyond which the 
excitation process is nearly uncorrelated. The time scale of the response is the measure of the memory 
duration of the system which is generally about one quarter of the natural period of a mode which contributes 
significantly to the total response. 

The Fokker-Planck equations for the rotor-bearing system are solved, to obtain the first order distribution of 
the response, which is further processed for the inverse problem of parameter estimation. The linear and 
nonlinear bearing stiffness parameters are obtained. Tire technique has a distinct advantage that, it does not 
require an estimate of the excitation forces and works directly on the response signals from the bearing caps. 
In this analysis the rotor is treated as a rigid body. The algorithm is tested by Monte Carlo simulation. The 
procedure is illustrated for a laboratory rotor rig and the results are compared with those from the analytical 
formulations of Harris (1984) and Ragulskis et al. (1974). 

The analysis becomes more involved, if the shaft flexibility is taken into account. The present study further 
attempts the problem of bearing stiffness parameter estimation for flexible rotors. In contrast to the rigid 
rotor case, which could be treated as a single degree freedom problem, the flexible rotor poses a nonlinear 
multi -degree of freedom problem. The excitation to the balanced system is taken to be random in nature, 
primarily arising out of bearing and assembly imperfections. The problem is formulated for a flexible shaft 
carrying a centrally located rigid disc and supported at the ends by nonlinear bearings. The inverse problem 
is approached by initially effecting a coordinate transformation, so as to enable the governing equations to be 
modeled as Markov Processes through the Fokker-Planck equations. The solution to the Fokker-Planck 
equation is obtained under certain engineering assumptions. A curve fitting algorithm is developed to process 
the statistical response of the system obtained by the solution of the Fokker-Planck equation, to extract the 
rotor-bearing stiffness parameters. 

The problem, for rotors with flexible shafts and carrying more than one disc, is considered next. The 
governing nonlinear differential equations for such multi-mass flexible systems are derived for balanced 
rotors with random excitation at the bearings. The governing equations are subjected to a coordinate 
transformation and modeled as Markov Processes. General form expressions, for the first order probability 
statistics of the response are obtained. The statistical response is processed to extract the rotor-bearing 
stiffness parameters. The procedure is illustrated for a laboratory test-rig with two discs and the 



experimental results are compared with the analytical guidelines of Harris (1984) and Ragulskis et al. 
(1974). The algorithm developed is tested by Monte Carlo numerical simulation procedure. 

The above studies concern single-degree-of-freedom and multi-degree-of-freedom systems, where the rotor is 
assumed to be balanced. The study, further, explores the possibility of estimation of linear and non-linear 
stiffness parameters, for rotors with unknown unbalance. The problem is formulated as single degree of 
freedom system. The case of a rotor with a rigid shaft, in nonlinear flexible bearings is investigated. The 
excitation to the system consists harmonic forces due to the unbalance and random forces due to arbitrary 
deviations, of bearing contact surfaces and subsurfaces, from their ideal design and their progressive 
deterioration during operation. These random forces are comparable to the harmonic excitation forces, if the 
unbalance is not significantly large. The parameter estimation procedure is based on the averaging technique 
of Bogoliubov and Mitropolsky (1961) for deterministic non-linear systems, extended by Stratonovich (1967) 
for stochastic differential equations. The governing equation of motion is transformed from the rapidly 
varying variables, namely displacement and velocity, to variables, amplitude and phase, varying slowly with 
time. Stochastic averaging is done, to take into account the effect of the random excitation multiplied by a 
correlated term, so as to model the slowly varying amplitude as an approximate Markovian process. A 
second order stochastic approximation is carried out and a one-dimensional Fokker-Planck equation is 
derived to describe the Markovian amplitude process. The response to the Fokker-Planck equation is derived 
and processed for parameter estimation. Along with the bearing stiffness parameters, estimates of the 
unknown unbalance of die rotor, its angular location and die damping ratio are also obtained, as by-products. 
The procedure is illustrated for a laboratory rotor rig and the experimental results are validated through 
comparisons with the available analytical guidelines. 

To summarize, non-linear bearing stiffness estimation procedure, based on statistical methods, have been 
developed for cases of rigid rotors, single disc flexible rotors, multi-disc flexible rotors. The study is 
extended, to include harmonic excitation to the non-linear system along with random excitation and the case 
of an unbalanced rigid rotor is discussed. The estimation procedure has an advantage over existing ones, for 
it does not require an estimate of the excitation forces and works directly on the measured response signals of 
the system and therefore can be employed on-line. The algorithms are illustrated for a laboratory rotor- 
bearing test rig and the results are compared with those obtained through an existing analytical model. 
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CHAPTER 1 


INTRODUCTION 


Bearing analysis constitutes a major area in rotordynamic studies. The response and stability of a rot® 
system are critically dependent on bearing characteristics. It has become evident in recent years, that ar 
important class of rotor bearing phenomenon cannot be studied without adequately accounting for the 
nonlinear forces produced by the bearings. Bearing nonlinearity assumes a greater role for high speed anc 
low weight rotor applications, where the vibration amplitudes tend to be relatively large. While, dynamic 
characterization of fluid film bearings has drawn considerable attention from researchers, the information 
available, on the nonlinear aspects of rolling element bearings under dynamic conditions, is relatively scarce. 

Rolling element bearings, and in particular ball bearings, despite their mechanical simplicity, are known to 
display highly nonlinear behaviour and present some very complex rotor problems. The approximate nature 
of the analytical determination of these bearing properties is responsible for some of the unreliability in the 
prediction of the response and stability of a rotor system. Procedures (Ragulskis, et al., 1974; Harris, 1984) 
are available for estimation of bearing stiffness under static loading conditions in isolation of shaft. The 
existing parameter estimation techniques for overall rotor-bearing systems (e.g. Muszynska and Bently, 
1990) involve, a controlled input excitation , that needs to be given to the bearings. 

The present study attempts for the development of a parameter estimation procedure, for the non-linear 
elastic parameters of bearings, that can be used on-line , and is based on the analysis of easily accessible 
vibration signals, picked up from the bearing caps. 

The forces setting the rotor-bearing system into vibrations comprise of both - harmonic forces, due to rotor 
unbalance and random forces due to bearing surface imperfections, caused by random deviations from their 
standard theoretical design and progressive surface and subsurface deterioration. In addition excitation can 
be contributed by from random inaccuracies in the rotor-bearing-housing assembly etc. 

Stochastical dynamics methods provide an opportunity, to attempt the inverse problem of parameter 
estimation, in such a rotor-bearing system, governed by a non-linear equation of motion and having both 



deterministic and random excitations. Under certain engineering assumptions, the random excitation can be 
taken to be an ideal white noise process and the response of a dynamic system (linear or non-linear) can be 
modeled as a diffusive Markov Process or an approximate Markov Process. 

The structure of a Markov Process is completely determined for all future times by the distribution at some 
initial time and by a transition probability density function, which satisfies the Fokker-Planck-Kolmogorov 
(FPK) equation. The FPK equation can be solved to obtain analytical expressions for the first order 
probability distribution of the stationary response. For the response to be modeled by a Markov Process, it is 
necessary that the excitations be approximated by ideal white noises. This restriction can be, in principle, 
removed at the price of increasing the complexity of the system. Linear filters are introduced between the 
ideal white-noise excitations and the system, to produce realistic excitation spectra for the nonlinear system. 
The ideal white-noise excitations, driving the filters, guarantee that the response of the extended system is a 
Markov Process (of higher order), for which an extended FPK equation can be written. For such an 
approximate Markovian modeling stochastic averaging is applied, whereby rapid fluctuations are averaged to 
provide simpler equations for slowly fluctuating quantities. The procedure is a nontrivial extension of the 
Krylov-Bogoliubov (1937) averaging method for deterministic excitations, since it involves accounting for 
the averaged effect of a random excitation multiplied by a correlated response. 

The problem is initially attempted for balanced rotors, by considering the random forces from the bearings 
as the primary source of excitation, which are generally large enough to cause measurable level of vibrations. 
The dynamics of the rotor-bearing system is modeled as a Markov Process and Fokker-Planck equations are 
formulated. The drift and diffusion coefficients in the FPK equation are derived from the nonlinear equations 
of motion of the dynamic system (Crandall, 1958). The approach to the solution of equation is greatly 
simplified if the overall random excitation to the system, from the variety of sources, is treated as ideal white 
noise. While many engineering applications are based on this idealization, insofar as the excitation itself is 
concerned, it turns out that the response obtained through such models are quite acceptable if the time scale 
of excitation is much smaller than the time scale of the response (Lin, 1973). The time scale for the 
excitation is the correlation time, roughly defined as the length of time separation beyond which the 
excitation process is nearly uncorrelated. The time scale of the response is the measure of the memory 
duration of die system which is generally about one quarter of the natural period of a mode which contributes 
significantly to the total response. 



The Fokker-Planck equations for the rotor-bearing system are solved, to obtain the first order distribution of 
the response, which is further processed for the inverse problem of parameter estimation. The linear and 
nonlinear bearing stiffness parameters are obtained. The technique has a distinct advantage that, it does not 
require an estimate of the excitation forces and works directly on the response signals from the bearing caps. 
In this analysis the rotor is treated as a rigid body. The algorithm is tested by Monte Carlo simulation. The 
procedure is illustrated for a laboratory rotor rig and the results are compared with those from the analytical 
formulations of Harris (1984) and Ragulskis et al. (1974). 

The analysis becomes more involved, if the shaft flexibility is taken into account. The present study further 
attempts the problem of bearing stiffness parameter estimation for flexible rotors. In contrast to the rigid 
rotor case, which could be treated as a single degree freedom problem, the flexible rotor poses a nonlinear 
multi-degree of freedom problem. The excitation to the balanced system is taken to be random in nature, 
primarily arising out of bearing and assembly imperfections. The problem is formulated for a flexible shaft 
carrying a centrally located rigid disc and supported at the ends by nonlinear bearings. The inverse problem 
is approached by initially effecting a coordinate transfonnation, so as to enable the governing equations to be 
modeled as Markov Processes through the Fokker-Planck equations. The solution to the Fokker-Planck 
equation is obtained under certain engineering assumptions. A curve fitting algorithm is developed to process 
the statistical response of the system obtained by the solution of the Fokker-Planck equation, to extract the 
rotor-bearing stiffness parameters. 

The problem, for rotors with flexible shafts and carrying more than one disc, is considered next. The 
governing nonlinear differential equations for such multi-mass flexible systems are derived for balanced 
rotors with random excitation at the bearings. The governing equations are subjected to a coordinate 
transformation and modeled as Markov Processes. General form expressions, for the first order probability 
statistics of the response are obtained. Tire statistical response is processed to extract the rotor-bearing 
stiffness parameters. The procedure is illustrated for a laboratory test-rig with two discs and the 
experimental results are compared with the analytical guidelines of Harris (1984) and Ragulskis et al. 
(1974). The algorithm developed is tested by Monte Carlo numerical simulation procedure. 

The study, further, explores the possibility of estimation of linear and non-linear stiffness parameters, for 
rotors with unknown unbalance. The problem is formulated as single degree of freedom system. The case of 
a rotor with a rigid shaft, in nonlinear flexible bearings is investigated. The excitation to the system consists 



harmonic forces due to the unbalance and random forces due to arbitrary deviations, of bearing contact 
surfaces and subsurfaces, from their ideal design and their progressive deterioration during operation. These 
random forces are comparable to the harmonic excitation forces, if the unbalance is not significantly large. 
The parameter estimation procedure is based on the averaging technique of Bogoliubov and Mitropolsky 
(1961) for deterministic non-linear systems, extended by Stratonovich (1967) for stochastic differential 
equations. The governing equation of motion is transformed from the rapidly varying variables, namely 
displacement and velocity, to variables, amplitude and phase, varying slowly with time. Stochastic averaging 
is done, to take into account the effect of the random excitation multiplied by a correlated term, so as to 
model the slowly varying amplitude as an approximate Markovian process. A second order stochastic 
approximation is carried out and a one-dimensional Fokker-Planck equation is derived to describe the 
Markovian amplitude process. The response to the Fokker-Planck equation is derived and processed for 
parameter estimation. Along with the bearing stiffness parameters, estimates of the unknown unbalance of 
the rotor, its angular location and the damping ratio are also obtained, as by-products. The procedure is 
illustrated for a laboratory rotor rig and the experimental results are validated through comparisons with the 
available analytical guidelines. 



6 


The drift and diffusion coefficients in the FPK equation can be derived from the nonlinear equations of 
motion of the dynamic system (Crandall, 1966). A general, closed form solution to FPK equation is yet to be 
found. When the time derivative in FPK equation is set equal to zero, the equation describes the first order 
probability distribution for the stationary response. The distribution is readily obtained for a limited class of 
vibratory systems with nonlinear restoring forces and special forms of nonlinear damping (Fuller, 1969; 
Caughey, 1971; Caughey and Ma, 1983; Pradlwarter et al., 1991). The Fokker-Planck equation was derived 
and the exact stationary response was obtained for certain cases of two-degree of freedom non-linear 
dynamic systems, by Ariaratnam (1960). The theory was generalized later by Caughey (1963) for multi- 
degree of freedom cases. 

The exact mean-square response and average, frequently provided by the first-order distribution, have been 
widely used as a basis for comparison with results of various approximate methods and/or procedures for 
computer simulation. The time dependent solution to FPK equation i. e., the transition probability density 
function, has been obtained for only a few nonlinear systems (Caughey, 1971; Roberts, 1981a). Yong and 
Lin (1987). Lin and Cai (1988a,b) have developed a systematic procedure to obtain the exact stationary 
response for either external (additive) or parametric (multiplicative) excitations, or both. Exact stationary 
solution in terms of probability density function, for a class of non-linear systems driven by a non-normal- 
delta-correlated process has been obtained recently by Vasta (1995). 

2.2 STOCHASTIC AVERAGING OR APPROXIMATE MARKOV PROCESS APPROACH 

For the response to be approximated by a Markov process, it is necessary that the excitation be 
approximated by ideal white noise. This restriction can be, in principle, removed at the price of increasing the 
complexity of the system, by introducing linear filters between the ideal white-noise excitations and the 
system. The filters serve to produce realistic excitation spectra for the nonlinear system, while the ideal 
white-noise excitations driving the filters guarantee that the response of the extended system is a Markov 
process (of higher order), for which an extended FPK equation can be written. This strategy has often been 
suggested but it has, apparently, never been implemented for nonlinear response problems. However, it can 
be shown by using stochastic averaging principles that, under certain conditions, the response of a nonlinear 
dynamic system to non-white excitation can be approximated by a Markov process. The designation - 
stochastic averaging - has been applied to a class of procedures in which rapid fluctuations are averaged to 
provide simpler equations for slowly fluctuating quantities. The procedure is a nontrivial extension of the 
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Kiylov-Bogoliubov (1937) averaging method for deterministic excitations, since it involves accounting for 
the averaged effect of a random excitation multiplied by a correlated response. 

There are three distinct approaches to stochastic averaging. The first approach was introduced by 
Stratonovitch (1967) and rigorously justified by Khasminskii (1966) and by Papanicolaou and Kohler 
(1974). Approximate equations for the slowly varying quantities are obtained by time-averaging the rapid 
fluctuations. This method, called the standard stochastic averaging method (or Stratonovitch's stochastic 
averaging method), is applied to narrow-band responses which can be represented as sinusoidal oscillations 
with slowly varying amplitude and phase. The formulae of the stochastic averaging methods were also 
derived by Lin (1986) using a different procedure, making the physical implication of the method clearer and 
more appealing to engineers. A generalization of stochastic averaging in random excitation with special 
emphasis was placed on casting the problem in a more formal mathematical framework by Red-Horse and 
Spanos (1992). Zhu et. al. (1994) improved the stochastic averaging procedure. Bouc (1994) has obtained 
the random response with a higher order approximation in the presence of a large non-linear stiffness term in 
a dynamic system subjected to random excitations. 

Another approach involves averaging the drift and diffusion coefficients in the FPK equation, with respect to 
time. The mathematical basis for this method is contained in a theorem by Khasminskii (1963) and is called 
the Averaging Method of Coefficients in FPK equation. One more approach is based on dividing the response 
variables into rapidly varying quantities and slowly varying quantities and approximate equations for the 
latter are obtained by averaging the rapid fluctuations of the former. The basis for this is also due to 
Khasminskii (1968). For a single degree of freedom oscillator, the rapidly varying quantity is the 
displacement and the slowly varying quantity is the energy envelope (Dimentberg, 1980; Roberts, 1983; 
Spanos, 1983; and Zhu, 1983). This method is called the Generalized Stochastic Averaging Method (or 
Stochastic Averaging Methods of Energy Envelope). Moshchuk et al. (1995) gave an analytical method 
based on the stochastic averaging of the energy envelope to treat the dynamic behaviour of single-degree-of- 
freedom elastic ocean structures subjected to non-linear hydrodynamic loading by Gaussian ocean waves. 

All the above three approaches have been extensively applied to nonlinear random vibration of mechanical 
and structural systems as effective approximations for prediction of the response, stability and reliability. 
Among the three stochastic averaging methods, the Standard Stochastic Averaging Method finds the most 
popular application on account of its accuracy, applicability and simplicity (Zhu, 1988). For nonwhite 
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excitation in non-linear random dynamic system, Cai (1995) developed a procedure of quasi-conservative 
averaging. He solved the Duffing oscillator under both external and parametric excitations with nonwhite 
spectral densities. 


2.3 STATISTICAL LINEARIZATION METHOD 

A natural method of approaching a nonlinear problem is to replace the given set of nonlinear equations by an 
equivalent set of linear ones, since a linear system is so much easier to analyse. This is achieved by 
minimizing the difference between the non-linear and linear sets of equations in an appropriate sense. The 
basic development of a suitable linearization procedure (variously known as statistical linearization, 
equivalent linearization or stochastic linearization), for randomly excited nonlinear systems, is usually 
attributed to Booton et al (1953), Booton (1954) and Caughey (1963a). They generalized the deterministic 
linearization methods of Krylov and Bogoliubov (1937) to the stochastic case. 


Although, rationalizations for the procedure usually appeal to a small nonlinearity argument, it fortuitously 
turns out that the error in mean-square response remains quite modest even for large nonlinearities 
(Jazwinski 1970). A great advantage of the statistical linearization method is that it can easily be generalized 
to cope with multi-degree of freedom systems, including those where hysteretic elements are incorporated. 
Tie earliest extensions of the theory in this direction were given by Caughey (1963a) and Kazakov ( 1 965a , 
1965b). Subsequently there have been a number of theoretical advantages in this area (Roberts 1981b). The 
technique has no difficulty in dealing with non-white excitations and can be further generalized to cope with 
non-stationary excitations and responses (Spanos 1981a, Roberts 1981b; Roberts and Spanos, 1990). Chang 
and Young (1989) obtained the stationary response of robot manipulators under stochastic base and external 
excitations by statistical linearization method. Socha and Soong (1991) reviewed and gave an assessment of 
the procedures in the area of statistical and equivalent linearization in the analysis of non-linear stochastic 
systems. Stochastic linearization on nonlinear multi-degree-of-freedom systems under random parametric 
excitation was examined by Falsone (1992). Elishakoff and Cai (1993) obtained an approximate solution for 
nonlinear random vibration problems by a partial stochastic linearization. Stochastic equivalent linearization 
was used by Chang (1994) to develop a finite element formulation for the dynamic response analysis of non- 
linear hysteretic plates subjected to random excitations. Soize (1995) described the calculation of the power 
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spectral density function of stationary response and the method of identification, of single degree of freedom 
non-linear second order dynamic systems excited by a white or a broad-band Gaussian noise. Grigoriu 
(1995) has developed an equivalent linearization method for the analysis of the non-linear systems with 
Poisson white noise input. 

2.4 PERTURBATION METHOD 

Perturbation methods are generally employed, when the amount of nonlinearity in a system is controlled by a 
small scaling parameter. The solution is sought in the terms of a power series in the small scaling parameter 
and successive terms are evaluated as linear responses to nonlinear functions of the preceding terms. This 
classical approach for detenninistic nonlinear problems (Stoker, 1950) was extended to random vibration 
problems by Crandall (1963). In practice, the calculations are seldom carried beyond the first perturbation. 
The perturbation series is often an asymptotic series, so that the higher order perturbations improve the 
approximation for small scaling parameter at the expense of worsening it for large scaling parameter. 

Functional series methods offer an alternative approach to developing an expansion, based on the linear 
solution. An example of the application of a such method is the work of Orabi and Ahmadi (1987). They 
used a Wiener-Hermite expansion and presented a formal procedure for deriving the deterministic equations 
governing the kernel functions, arising in the expansion (Roy and Spanos, 1990). A common difficulty with 
all expansion methods lies in establishing the regimes of convergence, in the appropriate parameter space. It 
is frequently found that, due to a combination of poor convergence properties and excessive computational 
requirements, it is only possible to obtain reliable results if there is a very small degree of nonlinearity. 

2.5 METHOD OF MOMENTS 

A set of differential equations for various statistical moments, or related quantities known as cumulants (or 
semi-invariants) and quasi-moments (Stratonovitch, 1967) of the response, as the function of time, can be 
obtained by multiplying the FPK equation by suitable functions and integrating over the probability space. 
Equivalent sets of equations can be derived directly from the dynamic equations of motions or the equivalent 
Ito equations. The response problem is thus reduced to a set of coupled ordinary differential equations. For 
linear systems the set of equations closes, in the sense that for some finite ‘n’ there is a set of V independent 
equations in which only n’ unknown moments appear. For nonlinear systems, the set of equations forms an 
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infinite hierarchy: There are always more unknown moments than equations. Approximate solutions have 
been proposed, based on ad-hoc closure assumption: Certain higher order moments are neglected or assumed 
to be related to lower order moments, in the same way that moments of Gaussian processes are related 
(Bolotin, 1979; Ibrahim and Roberts, 1978; Roberts, 1981a). 

A non-Gaussian closure method, based on constructing an approximate representation of the system’s 
probability density function with ‘n’ adjustable parameters and using ‘n’ independent moment equations to 
determine the parameters, was proposed by Dashevskii and Liptser (1967). A characteristic feature of this 
approach is that the complexity of the moment equations dramatically increases as the order of closures 
increases. However, results for nonlinear oscillators have been obtained for ‘n’ values up to 6, using both 
cumulant closure (Wu and Lin, 1984; Lin and Wu, 1984) and quasi-moment closure (Bover, 1978). These 
studies demonstrate that a significant improvement in accuracy can be obtained by progressing beyond 
simple Gaussian closure. The method has been applied by Ibrahim and coworkers ( Ibrahim, 1985; Ibrahim 
and Soundararajan, 1985; Ibrahim et al, 1985) to a study of the response of systems with nonlinear inertial 
terms to random excitation. The method gives good results in some cases but is apparently not a very robust 
procedure ( Crandall, 1985). 

A second-order closure method was presented by Nayfeh and Serhan (1990) for determining the response of 
non-linear systems to combined deterministic and random excitations. A closure method was proposed by 
Grigoriu (1991) for calculating moments of the state vector of non-linear system satisfying an Ito stochastic 
differential equation. Paola et al. (1992) presented the extension of the Ito rule for the case of vector real 
functions of the response of nonlinear systems excited by zero-mean Gaussian white noise processes. Davies 
and Liu (1992) have obtained the response of a non-linear oscillator excited by random narrow-band noise by 
stochastic averaging followed by partial linearization. Ibrahim et al (1993) have examined stochastic 
bifurcation in moments of a clamped-clamped beam response to wide band random excitations analytically 
(Gaussian and non-Gaussian closure), numerically (Monte Carlo simulation) and experimentally. An 
efficient moment calculation method has been developed, recently, by Katafygiotis and Beck (1995). 

2.6 METHOD OF EQUIVALENT NON LINEAR EQUATIONS 

An alternative generalization of statistical linearization has been proposed by Caughey (1986). The idea is to 
replace the original set of nonlinear differential equations by an equivalent nonlinear set, where the latter 
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belong to a class of problems which can be solved exactly. This class is, at present, very limited, and thus the 
range of applicability of the technique is correspondingly restricted. The results have been obtained for 
oscillators with nonlinearity in damping and stiffness (Caughey, 1986; Cai and Lin, 1988; Zhu and Yu, 
1989). The equivalent statistical quadratization methods have been applied to non-linear multi-degree-of- 
freedom systems subjected to random excitation by Spanos and Donley (1991, 1992). It has been 
demonstrated that the method is very effective as a means of predicting the probability distribution of the 
response, with reasonable accuracy. 

2.7 METHOD OF COMPUTER SIMULATION 

Numerical simulation or the Monte Carlo method (Shinozuka, 1972; Bolotin, 1979; Spanos, 1981b; 
Rubinstein, 1981; Spanos and Mignolet, 1989) consists of generating a large number ‘m’ of sample 
excitations, computing the corresponding response samples and processing them to obtain the desired 
response statistics. The backbone, of any digital simulation study, is an algorithm which provides a set of 
pseudo-random numbers, belonging to a population with a specified probability density function. Proper 
processing of this set of numbers can yield the values of sample functions of random process excitations, 
with pre-selected frequency content and temporal variation of intensity, at successive discrete equi-spaced 
times. Upon generating a single sample of excitation, the commonly available subroutines for the numerical 
integration of the differential equations can be employed to obtain the system response. Another sample of 
the excitation can then be generated and the computed values of the system response can be used to update 
its statistics. The procedure is, in principle, very general and applicable to stationary or non-stationary 
response of systems of any degree of complication. 

The statistical uncertainty in the response statistics decreases in proportion to ‘m’ while the computational 
cost increases essentially in proportion to £ m\ To gain one additional significant figure in a result requires a 
hundred fold increase in computational cost. Spanos (1981a) has estimated that, for cases where both 
statistical linearization and numerical simulation can be applied, the computational efficiency of the former 
will be of the order of 100 to 1000 times better than the latter. In applications where no exact solution is 
available numerical simulation is often used as a basis for assessing the accuracy of other approximate 
methods. Benaroya and Rehak (1988) reviewed the field of probabilistic structural analysis where finite 
element methods are used. The stochastic finite element analysis was applied to obtain response of 
geometrically and materially non-linear plane trasses under random excitations by Chemg and Wen (1991). 
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The finite element method is applied, in conjunction with the method of equivalent linearization, to large - 
deflection random response of thermally buckled beams by Locke and Mei (1990) and to the non-linear 
random response of beams by Locke (1994). 

2.8 ROLLING ELEMENT BEARING VIBRATIONS AND STIFFNESS ESTIMATION 

Estimation of the elastic parameters of bearings involves establishing a relationship between the incident load 
on the bearing and its resultant deformation. The classical solution for the local stress and deformation of 
two elastic bodies apparently contacting at a single point was established by Hertz in (1896). Hertz's analysis 
is applied to surface stresses caused by a concentrated force, applied perpendicular to the surface. In the 
determination of contact deformation versus load, the concentrated load applied normal to the surface alone, 
is considered, for most rolling element bearing applications. Methods of calculation of the surface and 
subsurface stresses of the combination of normal and tangential (traction) stresses are complex, (Zwirlein 
and Schlicht, 1980). Owing to infinitesimally small irregularities in the basic surface geometries of the 
rolling contact bodies, neither uniform normal stress field nor a uniform shear field are likely to occur in 
practice (Sayler et al, 1981; Kalker, 1982). Rigorous mathematical/numerical methods have been developed 
to calculate the distribution and magnitude of surface stresses in any line contact situation, that is, including 
the effects of crowning of rollers, raceways, and combinations thereof (Kunert, 1961; Reusner, 1977). 
Additionally, finite element methods (FEM) have been employed (Fredriksson, 1980) to perform the same 
analysis. 

It is possible to determine how the beanng load is distributed among the balls or rollers, after having 
determined how each ball or roller in a bearing carries load. To do this it is necessary to develop load- 
deflection relationships for rolling elements contacting raceways. Most rolling bearing applications involve 
steady-state rotation of either the inner or outer raceways or both. Rolling element centrifugal forces, 
gyroscopic moments and frictional forces and moments do not significantly influence this load distribution in 
most applications. Theoretical models (Palmgren, 1959; Ragulskis et al., 1974; Harris, 1984; Eschamann et 
al. 1985; Stolarski, 1990) are available for estimation of bearing stiffnesses under static loading conditions. 

Bearings concern vibrations caused due to geometric imperfections of contact surfaces were first analyzed by 
Lohman (1953) and Gustavsson (1962). Significant contributions have been made by Gupta et al (1975, 
1977, 1978a, 1978b, 1978c, 1978d) towards the understanding of rolling element bearing dynamics. 
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Comprehensive investigations have bear carried out on the high frequency response of bearings ( McFadden 
and Smith, 1984) and its relation to surface irregularities ( Sunnesjo, 1985; McFadden and Smith, 1985; Su 
et al, 1993). Lim and Singh have analyzed the vibration transmission through rolling element bearings in 
series of publications (1990a, 1990b, 1991, 1992, 1994). 

A method for determination of the nonlinear characteristics of bearings using the procedure of Krylov- 
Bogoliubov-Mitropolsky has been suggested by Kononenko and Plakhtienko ( 1970). Honrath (1960) and 
Elsermans et al (1975) examined the stiffness and damping of rolling element bearing experimentally. 
Walford and Stone (1980) designed and fabricated a test rig for direct measurement of the relative 
displacement of the shaft and bearing housing for the oscillating force applied to the bearing housing, which 
is used to obtain the stiffness parameters. They found that the interfaces play between races and housing and 
shaft play a significant role in the determination of bearing stiffness and damping. Beatty and Rowan (1982) 
determined the dynamic stiffness of ball bearing. Kraus et al (1987) presented a method for the extraction of 
rolling element bearing stiffness and damping under operating conditions. The method is based on 
experimental modal analysis combined with a mathematical model of the rotor-bearing-support system. 
Effects of speed, preload and free outer race bearings on stiffness and damping have been investigated. The 
technique involves a controlled input excitation to be given to the bearings. Mitchell et. al (1966) obtained 
the stiffness of an oil film bearing, experimentally, by application of static loads. Morton (1971) devised the 
measurement procedure for estimation of the dynamic characteristics of a large sleeve bearing by application 
of dynamic loads using vibrators. Nordmann and Schollhom (1980) identified the stiffness and damping 
coefficients of journal bearings by means of the impact method. Sahinkaya and Burrows (1984) estimated 
the linearised oil film parameters from the out-of-balance response where the shaft was excited by a known 
unbalance force. Muszynska (1990) has developed a perturbation technique for estimation of these 
parameters. Goodwin (1991) reviewed the experimental approaches to rotor support impedance 
measurement, with particular emphasis on fluid film bearing impedance measurement. A general procedure 
for identification of restoring force nonlinearity from a system’s response to a white-noise excitation has been 
discussed by Dimentberg and Sokolov (1991). Nonlinear stochastic contact vibrations and friction at a 
Hertzian contact have been studied by Hess et al (1992). Analytical and experimental studies are carried out 
by them, using the Fokker-Planck equation and simulating the vibrations to the contact region, either 
externally by a white Gaussian random normal load or internally by a rough surface input. Childs and Hale 
(1993) devised a test apparatus and facility to identify the rotordynamic coefficients of high speed 
hydrostatic bearings. 



CHAPTER 3 


BEARING STIFFNESS ESTIMATION IN SINGLE-DISC RIGID ROTORS 


The forces setting the rotor-bearing system into vibrations, comprise of both - harmonic forces, due to rotor 
unbalance and random forces due to bearing surface imperfections, caused by random deviations from their 
standard theoretical design and progressive surface and subsurface deterioration. The problem is initially 
attempted for balanced rotors, by considering the random forces from the bearings as the primary source of 
excitation, which are generally large ©rough to cause measurable level of vibrations. 

3.1 EQUATION OF MOTION 

The governing equation for a balanced rigid rotor supported at ends in bearings (Figure 3.1) with nonlinear 
stiffnesses can be written as 

x+2gco n x+60n\_x + XG(x)] = f(t) (3.1) 

where G can be a polynomial in x and X is the unknown non-linear stiffness contribution parameter. f{t) 
in equation (3.1) represents the random excitation to the system caused due to bearing surface 
imperfections and progressive wear and tear. In addition excitation can be contributed by from random 
inaccuracies in the rotor-bearing-housing assembly etc. 

The Markov vector approach is adopted for obtaining expressions for the response of the nonlinear random 
system of equation (3.1). The approach to the solution of equation (3.1) is greatly simplified if the overall 
random excitation to the system, from the variety of sources, is treated as ideal white noise. For an ideal 
white noise process, the response of a dynamic system (linear or non-linear) is a diffusive Markov process. 
The structure of a Markov process is completely determined for all future times by the probability 
distribution at some initial time and by a transition probability distribution, which satisfies the Fokker- 
Planck equation. While many engineering applications are based on this idealization, insofar as the 
excitation itself is concerned, it turns out that the response obtained through such models are quite 




Figure 3.1 


Rotor on rolling element bearings 
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acceptable if the time scale of excitation is much smaller than the time scale of the response (Lin, 1973). 
The time scale for the excitation is the correlation time, roughly defined as the length of time separation 
beyond which the excitation process is nearly uncorrelated. The time scale of the response is the measure of 
the memory duration of the system which is generally about one quarter of the natural period of a mode 
-which contributes significantly to the total response. The statistical averages of the ideal white noise can be 
expressed as 


</( ti )>=0 

</('i)/(' 2 )>= 2 nS 0 5(t 2 -t x ) 


(3.2) 


where S 


is the uniform spectral density of excitation and 8 denotes Direc delta function. 


3.2 F-P-K EQUATION 


For a set of state space equations 


ex, / ct = ir,. 

&, / a = f(t ) - (1 / Mj - 8^( Xl ,x 2 


i = 1> 2, ... JV 


(3.3) 


with fi(t), being the zero mean white noise excitation, as defined in equations (3.2), the response is a 
diffusive Markov process, for which the transition probability density function is governed by the Fokker- 
Planck equation, given as (Caughey, 1963) 


cpi dt = ~Y\dl cx t {a : p} + 8/ ox. {b,p}] 

/=l 

+ (l/2)XZt^ 2 ;&f{c ii p}+d 2 / cx] {d JJ p} + d 1 !& t &j{e v p}] 

i=i j=i 


(3.4) 


where N is the number of degrees of freedom in the dynamic-system-model. 
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The first order moments, a j , b { , in the Fokker-Planck equation, called the drift coefficients, are given by 
a, (x, t) = lim (< Ax, > /At) 

Ar-*-0 

/ = 1,2, ... ,N (3.5) 

bAx, r) = lim (<Ax, >/At) 

A/-»0 1 


and the second order moments, c u , d M , e y , called the diffusion coefficients, are expressed as 
c„(x, t) = lim (< Ax, 2 > /At) 

A/-+0 

djj(x,t) = lim (< Ax 2 > /At) i = 1, 2, ... (3.6) 

e„ (x, x, t) = lim (< Ax, Ax , > /At) 

J AJ->0 J 


The governing equation, (3.1), for the rotor-bearing system under consideration can be rewritten as 


dx / cit = x 

dx / dt = fit)- 2C,co„x - af n x - <y 2 AG(x) 


(3-7) 


and its drift and diffusion coefficients, from equations (3.5) and (3.6), obtained as 


a, = x 

b x = -2£a> „x - <o 2 n x - co 2 n AG(x) 

and 

c u = 0 

d \ i = 2 ttS 0 

e n = 0 

The Fokker-Planck equation, for the rotor-bearing system, can then be written as 

dp d ? d“ p cb 

-x — +—( 2 S<»nXP) + a, n\. x + : L G ( x )'\-z; + 7:S o — T = IT 

dx. dx dx dx" dt 


(3.8) 


(3.9) 


(3.10) 
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3.3 RESPONSE 


For a stationary case, equation (3.10) reduces to 


dp d 2 dp d^ p 

— (2go) n xp) + (O n [x + AG(x)] -^ + nS Q — =- = 0 
Sc S Sc u 


for which the solution is 


p(x,x) = cexp x 2 +^-x 2 co 2 n +Ao) : „g(x)} ] 

flu Q 4 2t 


(3-11) 


where 

g{x) = ]a^)d^ (3.12) 

o 

The probability density functions p{x) and p(x) are obtained from equation (3.12) as (Roberts and 
Spanos, 1 990), as 


OO 

P(x)= \p(x,x)dx 

— OO 

= Cj exp [-^r-{l-x 2 +2g{x)}] 

2 

where 

c, = 7 / J exp [---^-{-x 2 + Xg(x))\dx 

-oo fl& 0 — 


and 


(3.13) 


P(x)= \p{x,x)dx 


1 \cco„ r 2cco n .1 , 2 s, 

«p[— VM t* }] 


it 


*S n 


(3.14) 



19 


The variance of the velocity response is 


oo 

CT- = \x : p(x)dx 


xSp 


(3.15) 


From equations (3.13) and (3.15), the probability density function for the displacement response can be 
written as 


co 


p(x) = c, exp [ H-x* +%(*)}] 

<Ji 2 


with 

oo 

c,=l/\ 


2 i 

exp [-^f {-x* + Xg{x)}]ctc 

<j\ 2 


(316) 


3.4 EXTRACTION OF BEARING STIFFNESS PARAMETERS 

The bearing stiffness parameters are obtained from the experimentally obtained random response in terms 
of the linear parameter a> 2 n and the nonlinear stiffness parameter X . These parameters are obtained for 
both, the vertical and horizontal directions, the problem formulation, in the horizontal direction, remaining 
identical to that in the vertical direction. (The stiffness cross-coupling in rolling element bearings being 
negligible, the governing equations in the vertical and horizontal directions are taken to be independent, 
during the entire course of the present study.) 

Tire probability density function for any two displacements x, and x M (x (+; > x,), from equation (3.16) 
are 

2 i 

P(Xj) = C/ exp [-— t{~ 2 x \ +4g(*,)}] 

2 


(3.17) 
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and 


P( x j + ]) = c l expl- 


oit 


N 2 


i+1 


+ %(*,+/)}] 


Defining 
A*,- = x m - x - 

so that 

*1+1 = (Y, + AY,) 2 
» x? + 2x i Ax j 
for small Ax. 


and 

gix i+l ) = g(x i +Ay,.) 


a g(Xi) + ^^-[ x=Xi Ay, 


and substituting from equations (3.20) and (3.21) into equation (3.18), gives 


P(x m ) = c, exp [- %• x! + Xg(x . )] exp [- %{x, + X }Ax, ] 

CT ■ ^ O" i C3DC 


Combining equations (3.22) and (3.17) gives 

P(*m) = Pi*i)exp [-%{*, + } Ax J 

a, ax 


(3.18) 


(3.19) 


(3.20) 


(3.21) 


(3.22) 


(3.23) 


For N displacement values, x ; ,x 2 , ,x N , equation (3.23) can be expressed as a set of N-l linear 

simultaneous algebraic equations, as. 


A*, P( X m) <°n 
i = l,2 N-l 


dgjx) , 

dx lx=x ' 


X = x, 


(3.24) 
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Equations (3 .24) are solved for 1 /co 2 and X , using the least square fit technique. The variance, <j\ and 
the probability function, p(x), are computed from the experimentally obtained displacement and velocity 
data ( x and x ), which are taken as zero mean Gaussian processes and the nonlinear spring force provided 
by the rolling element bearings is taken to be cubic in nature i.e. G(x) = x 3 , based on the work of 
Ragulskis et al (1974). They conducted exhaustive experimental work on isolated ball bearing and found 
that cubic nonlinearity best fits their experimental data. Subsequently several other researchers have also 
employed the cubic nonlinearity model. 

3.5 EXPERIMENTATION 

The laboratory rig for the experimental illustration of the technique is shown in Figures 3.2, 3.3 and 3.4. 
The rig consists of a disc centrally mounted on a shaft supported in two identical ball bearings. The shaft is 
driven through a flexible coupling by a motor and the vibration signals are picked up (after balancing the 
rotor, balancing is done at 1000 rpm and vibration signals are picked up at that speed only) in both, the 
vertical and horizontal directions, by accelerometers mounted on one of the bearing housings. The signals 
from the accelerometers are digitized on a PC/AT (Figure 3.5) after magnification. (The specifications of 
the instrumentation is given in the Appendix - A) 

Typical displacement and velocity signals, in the vertical direction, picked up by the accelerometer are 
given in Figures 3.6 and 3.7. The probability density function, p(x), of the displacement is shown in Figure 
3.8. Tire bearing parameters estimated from equations (3.24) are - 

Vertical direction: 

a> l = 5.42 x 10 ' (rads/sec)' 

X = - 1.27 x 10* (mm -2 ) 

Horizontal direction: 

co l = 3.21 x 10 ' (rads/sec)' 

X =- 1.29 x 10 5 (mm -2 ) 

(Refer to Appendix C for statistical Error Table) 
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Figure 3.5 


Digital storage oscilloscope connected to PC/AT with GPLB 



Velocity x (mm/sec X 10 1 ) Displ. x (mm X 10 4 ) 



Figure 3.6 Displacement signal in the vertical direction 



Time t (sec X 10 1 ) 

Figure 3.7 Velocity signal in the vertical direction 



Figure 3.8 Probability density distribution 

of vertical displacement 





3.6 MONTE CARLO SIMULATION 


The algorithm is tested by Monte Carlo simulation, for the assumption involved. The values of oy n and X , 
estimated from the experimental data, are fed into equation (3.1). A broad band excitation force, f(t), with 
zero mean and Gaussian probability distribution, typically described in Figures 3.9 and 3.10 is 
computationally simulated. This force is also fed into equation (3.1) and the displacement and velocity 
responses, x and x , are numerically obtained through fourth order Runga-Kutta numerical technique. This 
simulated vertical displacement and velocity response is shown in Figures 3.11 and 3.12. The probability 
distribution of the simulated vertical displacement is shown in Figure 3.13. Parameter estimation, as 
described in the above subsection is now carried out with the simulated response, by feeding it into equation 
(3.24) to obtain the values of co~ n and X . A similar exercise is carried out to obtain the parameters in the 
horizontal direction. These values are listed in Table 3.1. 

Table 3.1 Estimated and Simulated Bearing Stiffness Parameters 



Parameters 


Experimental 

Simulated 


<ZT (rads/sec)' 

X (mm ' ) 

C0 2 n (rads/sec) 2 

X (mm 2 ) 

Vertical 

5.42 x 10 7 

- 1.27 X 10 15 

5.42 x 10" 

- 1.29 x 10* 

Horizontal 

3.21 x 10" 

- 1.29 x 10 6 

3.21 x 10" 

- 1.23 x 10* 


The good agreement between the values of the bearing stiffness parameters, o 2 n and X, obtained by 
processing the experimental data and those from the Monte Carlo simulation, indicate the correctness of the 
experimental and algebraic exercises. It should be noted that the simulated values of the bearing stiffness 
parameters are obtained for an ideal white noise excitation, while the experimental ones are obtained by 
processing the actual response of the system, where the unknown excitation was idealized as white noise. It 
also needs to be pointed out that the value of the damping ratio, g , is not required for the estimation 
procedure ( equation (3.24) ) . Any convenient value of g can be employed in equation (3.1) for the 
purpose of simulation ( g =0.02 has been assumed in the present simulation ). 




Probab. funct. p(f) Force / (N x 10 








Displ. x (mm X 10 4 ) 


Figure 3.13 


Probability density distribution 

simulated mtkal diaalattiasat 



19 


3.7 VALIDATION 

The values of the bearing stiffness parameters <X) : n and A , obtained by the procedure outlined are 
compared with those obtained from the analytical formulations of Harris (1984) and Ragulskis et al. 
(1974), which are based on Hertzian contact theory. 

For a given external radial load R on a bearing (Figure 3.14), the total elastic force at the points of contact 
of the rth ball with the inner and outer races is expressed as (Ragulskis et al., 1974) 

Si = K n (g + X cos Vi+y sin T ],) 2 ' 2 (3.25) 

and its projection along the line of action of the applied force is 

Ft = K^ig + xcosrjt + y sin Tj t) 3 ' 2 cost] t (3.26) 

where g is the radial preload or preclearance between the ball and the races and x and y are the 
displacements of the moving ring in the direction of the radial load and perpendicular to the direction of the 
radial load respectively, rj. is the angle between the lines of action of the radial load (direction of 
displacement of the moving ring) and the radius passing through the center of the rth ball. K n is a 
coefficient of proportionality depending on the geometric and material properties of the bearing. The 
specifications of the test bearing are given below. 


Ball beanng type 

SKF 6200 

Number of balls 

6 

Ball diameter 

6 mm 

Bore diameter 

10 mm 

Outer diameter 

30 mm 

Pitch diameter 

20 mm 

Inner ring ball race radius 

3.09 mm 

Outer ring ball race radius 

3.09 mm 

Allowable pre-load 

0 - 2 pm 

Rotor mass per bearing 

0.41 Kgs 
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The value of K n , for the test bearing with above specifications, is estimated by the method suggested by 
Harris (1984) as 2.82 xlO 5 


The total elastic force in the direction of the applied force is 

F = jrF t (3.27) 

1=7 

where n is the total number of balls in the bearing. 


Using the condition of zero elastic force in the direction perpendicular to the elastic load, the deformation, 
y, perpendicular to the radial force line is expressed as 


3 ; = 


^ [g + xcos( rjj )] 3/2 sin (77,) 

i=I 

n 

2 [g + xcos ( 7, )]" 2 sin 3 ( t], ) 

7=7 


(3.28) 


Equations (3.26) and (3.28) are used in equation (3.27) and the bearing stiffiiess is determined as a function 
of the deformation x as 

k(x) = dFj & (3-29) 

Substituting equation (3.27) in equation (3.29), taking into account equation (3.28) the bearing stiffness is 
expressed as a function of deformation as 


k{x) = K„ £ [g+x cos t]i - (A / Bn) sin 77, ] m [cos 77 , - {CBn - AD(n - 1)} / {Bn) 2 sin 77 ; ] cos 77, 
1 = 1 


where 

n 

^ = X[£ + xcos 7,] 3/2 sin 77 ,.; 

1=1 

n 

C = £[g + xcos77j 1/2 sin 77, cos ?7 f ; 

1=1 


B = 'L[g + xcosr li] sin 77 ,. 


7=1 


(3.30) 


D = Y,[g + XCOST 1 7 ] sin 2 77 ,. cos 77 , ; 

7=1 
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It can be seen that the bearing stiffness is critically dependent on the preloading, g, of the balls. While the 
manufacturer, may, at times, provide the preload range, the exact value of the preloading of the bearing 
balls in tire shaft-casing assembly, especially during operations which have involved wear and tear, would 
be difficult to determine. Tire stiffness of tire test bearing is plotted in Figure 3. 1 5 as a function of the radial 
deformation, x, for various allowable preload values, g. The bearing stiffness obtained experimentally, 
using the procedure developed, also shown in Figure 3.15, shows good resemblance to theoretically 
possible values. It is also to be noted that the theoretical stiffness calculations are based on formulations 
which analyse the bearing in isolation of the shaft. The comparison between the experimental and 
theoretically possible stiffnesses is also listed in Table 3.2. The expressions for the theoretical stiffnesses in 
Table 3.2 have bean obtained by curve fitting the stiffness values obtained from equation (3.30), through a 
quadratic in x. 


Table 3.2 Estimated and Theoretical (Ragulskis et al., 1974; Harris, 1984) 
Bearing Stiffness Parameters 


Preload 

(mm) 

Theoretical Stiffness 

(Radial) 

(N/mm) 

Estimated Stiffness 

(N/mm) 

0.0002 

1.20 x 10 '- 4.01 x 10 'V 




1.32 x 10 ^ - 5.08 x 10 10 x 2 (horizontal) 

0.0003 

1.47 x 10 ^ - 2.18 x 10 /0 x 2 





0.0004 

1.69 x 10'- 1.42 x 10 ; V 





0.0005 

1.89 x 10'- 1.02 x 10 7 V 




2.23 x 10'- 8.50 x 10 'V (vertical) 

0.0006 

2.08 x 10 ' - 6.09 x 10 9 x 2 



















Stiffness k(x ) (N/mm X 10 -4 ) 



Figure 3.15 Stiffness comparison 

1,2 - Present study (In vertical and horizontal directions respectively) 
3,4, 5,6,7 - Harris (1 984) and Ragulskis et al. (1974) with prelaod 
0.0002, 0.0003, 0.0004, 0.0005 and 0.0006 nun. respectively (In any 
radial direction) 




3.8 REMARKS 


The stiffness parameter estimates, from the procedure developed, show good agreement with the analytical 
formulations available for isolated ball bearings. The advantage of the proposed methods over other 
available techniques is distinct from the fact, that it does not involve measurement of the excitation forces 
and works directly on the random response signals, which can be conveniently picked up at the bearing 
caps and also that the procedure does not require the damping in the system to be known. 



CHAPTER 4 


BEARING STIFFNESS ESTIMATION 
IN FLEXIBLE SINGLE-DISC ROTORS 


The bearing stiffness estimation procedure, described in Chapter 3, becomes more involved, if the 
shaft flexibility is taken into account. In contrast to the rigid rotor case, which could be treated as a 
single degree freedom problem, the flexible rotor poses a nonlinear multi-degree of freedom 
problem. The excitation to the balanced system is taken to be random in nature, primarily arising 
out of bearing and assembly imperfections. The problem is formulated for a flexible shaft carrying 
a centrally located rigid disc and supported at the aids by nonlinear bearings 


4.1 EQUATIONS OF MOTION 


A balanced rotor, with a centrally located disc on a massless flexible shaft supported in bearings at 
ends is shown in Figure 4.1. The shaft is treated as free-free body, carrying unknown effective 
bearing masses m x and m 2 at its ends and the known disc mass m 3 in the center. The bearings are 
incorporated through external “forces”, F b , acting on masses m x and m 2 . Taking the shaft 
stiffness and damping forces as F s and F d , respectively, the equations of motion are written as 


-F di ~ F s, ~F„ x =m l x l 
~F d , -Fs, ~Fu = m 2 x 2 


(4.1) 

(4.2) 


-F, -F, =m,x. 


(4.3) 


The shaft stiffness forces are 
F s , = (^U X l ~^^12 X 2 "h ^13 X 3^ 

F Sl = (ft 12 X 1 + ^22 X 2 4~^-23 X 3 ) 

F s> = (^i3 x i 4"^23 X 2 4-A 33 x 3 ) (A4) 

The shaft stiffness parameter k tJ is defined as the rth force corresponding to a unity deflection with 
all other deflections held to zero and can be obtained from Strength of Materials formulae. 




Figure 4.1 Flexible rotor on rolling element bearings 


Bearing center line 
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Hie shaft damping forces are 
Fj r = («„*! +a n x 2 +a 13 i ;( ) 

Fj 2 — (^12-^1 + #22*2 + #23*3 ) 

F di = (# 13 *1 + #23*2 +«33 i 3) ( 4 ' 5 ) 

Shaft damping parameter, a {j , can be defined in a manner similar to k tj . 

The excitation to the system is taken to be random in nature, as in the previous chapter, mainly 
arising out of the bearing surface imperfections, caused by the random deviations from their 
standard theoretical design and progressive surface and subsurface deteriorations. Sources such as 
inaccuracies in the rotor-bearing-housing assembly etc. can also contribute to the random 
excitation. 

Taking 5, and s 2 as the effective random displacements at the bearings, primarily due to surface 
imperfections and inaccuracies in the rotor-bearing-housing assembly, the bearing forces on masses 
m x and m 2 can be written as 

F bl — {^r,, (*i — s i ) + k Ni t G(x, — ■?, )] 

F bl = {^(*2 ~ s z) + k- NL i G{x 2 - (4-6) 

In the above . k L and k NL are the unknown linear and nonlinear bearing stiffness parameters and 
G can be a polynomial in x . 

Equations (4.6) can be rewritten, more generally, as 

F b{ = { k L, + k NLl G ( x i )} “ F i (0 
F b2 — |^i 2 x 2 f k NLi G(x 2 )} — F 2 (t) 

where F x (t) and F 2 (t) are the random components of the bearing forces. 


(4.7) 
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Using (4.4), (4.5) and (4.7) the equations of motion (4. 1-4.3), can be written as 


r_ 

O 

O 

sf 

r .. 1 
x l 


i 

JP 

k> 

£> 

w 

1 


V 


dV {x x ,x 2J x 3 ) ! 3 c x 


'K 

0 m 2 0 

* X 2 

>4- 

a U a 22 a 2l 


x 2 

►4« 

oV ( ? x 2 yXj ) / dx*i | 

> = < 

?2> 

0 0 m 3 

4. 


13 ^23 


x 3 \ 


5V (X-^ , X 2 ? X 3 ) / ^3 j 


0 


where 

V(x t ,x 2 ,x 3 ) = m 3 

+nu 


+ l C ° 2 22 X l +C ° 2 \ 2 X 1 X 2 + ®ll X 2 X 2 + ^n X l X 3 

Z* jL* 


^< x ? +co\X l g{x l ) + ^co 2 n xl + <v 2 „A 2 g(x 2 ) 


(4.9) 


(4.10) 

(4.11) 

4.2 ORTHONORMAL TRANSFORMATION 

The Markov vector approach extended to nonlinear multi-degree-of-freedom systems (Nigam, 
1983) is adopted for the solution of equations (4.8). 

The equations of motion, ( 4 . 8 ), with damping and forces F x and F 2 set to be zero, are solved for 
eigenvalues p 2 , p 2 , p 2 and orthonormal modal matrix [£/] , so that 


*<*) = fcufw 


CD 


2 

m. 


, kr 

_,2 b 

= 


A, = 


1 NLj 


Q ) 2 nt, 
3 




38 


*n *11 * 11 1 \p] 0 0 

m T *,1 K * 1 . [£*]- 0 pi 0 

Jc n k 23 k 3 3 _ 0 0 p 2 _ 

m x 0 0 10 0 

[U] T 0 m 2 0 [£/]= 0 1 0 

0 0 m 3 0 0 1 


(4.12) 


(4.13) 


(The eigenvalue and eigenvector matrix elements are given in Appendix -B.) 


Application of coordinate transformation 



(4.14) 


to the equations of motion, (4.8), yields, 

T7-^(»7i,»7a ,*7 3 )/^7 i 

f .. -N P - I /' . 'I A/l ( 

1 0 0 Tjy /?„ 0 o V\ I \ ( h 

o l o i) 2 U o p 22 o \ij 2 \ = \q2 * 

.o o lJ^J [ o o 7 M 

— W(ti 1 ,ti 2 , tj 3 )/ drj 3 

l M 3 J 

where 

/?U 0 0 (X u CC 12 <Xj3 

0 p 22 0 =[(/] r a 12 a 22 a 23 [^] 

. 0 0 >^33 J L«13 a 23 «33. 


(4.15) 


(4.16) 
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^(7l,72>7 3 ) = ™3 


+ «3l7 2 +W2l7 3 ) 2 +^©L(»7 i +“32^ +«22 ? 7 3 ) 2 


1 j , 2 

+ 2 <y 33 ( ' 7 * +? ?2 +? ?3 ) 


+<y u(7i + “ 3 i 7 2 + u 2i7 3 )(7i + “ 3 27 2 +“ 22 ^ 3 ) 
+“>23 (7l +“32 72 +«22 7 3 )(7l +72 +7 3 ) 
+<y i 3 (7i +“ 31 72 +«2i7 3 )(7i +72 +73)] 


+m, 


-® 2 ,(7l +«3l72 +«2l7 3 ) 2 +* }2 , ;l lg(7l +«3l72 + «2l7 3 ) 


+«327 2 + “22 7 3 ) 2 +^ 2 2 X 2 g(r7 , +m 32 7 2 + « 22 7 3 ) 


f *s 

7 , 


X 

'72 



73 , 

■ 

0 


(4.17) 


(4.18) 


and the modal masses M, , M 2 and M 3 are given by 
M x = m x + m 2 + m 2 
M 2 = u\ x m x +m 2 + u] 2 m 3 

M 3 = u\ x m x +m 2 +u 22 m 3 (4.19) 

The approach to obtaining the response of the system is simplified, as in Chapter 3, by making an 
engineering assumption that random excitation to the system is such that the generalized forces q t , 
in equations (4.15) can be treated as ideal white noise. An argument, similar to that in Chapter 3, 
can be given regarding this idealization, that, it turns out that the response obtained through such 
models are quite acceptable if the time scale of the excitation is much smaller than the time scale of 
the response (Lin, 1973). The excitation of equation (4.15), is treated as uncorrelated Gaussian, 
white random forces with the following properties 


E[q x (t)] = 0; % 2 (0] = 0; % 3 (0] = 0 
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and 

-ti) 

^[^2 C^i )^ 2 C^ 2 )] = 2 /r^ 2 < 5 Ct 2 -t,) ( 4 . 20 ) 

^[93^1)93^2)] = 2 ^j 5 (t a -tj) 

where <t> x ,(f> 2 and ^ 3 denote the excitation intensity factors. 

4.3 F-P-K EQUATION 
Rewriting equations (4.15) as 


Vi +P„r] + (\I M i )^{ri l ,ri 2 ,r] i )l = o 


1 = 1,2, 3 (4.21) 


and then, in state space form as 


drj i I dt = t] { 

&Vi / dt = q,~ P n Vi -0/ M,. \dVtdv,) 


(4.22) 


the drift and diffusion coefficients (refer equations (3.5 - 3.7) ), for equations (4.22) can be written 
as 


a x = v x a 2 = i] 2 « 3 = V 3 

b x = <7, -PnVi -(1/ M l )[3'0h,V2,V 3 )S &h] 
b 2 =q 2 - p 22 V 2 - ( 1 !M 2 )W(v i,Vn V 3 ) 1 c?1 h ] 
b 3 = q 3 ~ Pa Vi ~ 0 / M 3 )[W(V 1 , V 2 , Vi ) / #Vi 1 


(4.23) 


and 

c u =0 c 22 =0 c n = ® 

d n = 27 t(j) x d 22 = ln<t> 2 d 33 = 2 n<j) 3 

e n ~ e 22 ~ e ii ~ e n = C 23 = e 31 = ® 


(4.24) 
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With the help of the above, the joint probability density function, p(J] l ,rj 2 ,Tj 3 ,rj l ,Tj 2 ,rf 3 ) , for 
the motion governed by equations (4.15) can be described by the Fokker-Planck equation 


-Vi 


dp x M x dp x dp x dp x 


1 dV dp d 


P u i] x p + 7t(f> x —- 

v dp u 


~V 2 


eft 1 dV eft 


dr] 2 M 2 dp 2 dp 2 


+- 


d_ 
dr] 2 


PriVi P + x</>7 


dp 




dp 


1 dV dp d 


dr] 2 ' dp 3 M 3 dp 3 dp 3 dp 3 


fi^riiP + x4> 3 


4? 


dpj 


eft 

dt 


(4-25) 


4.4 RESPONSE 


For a stationary case equation (4.25) reduces to 


. dp 1 dV dp d_ 

?h dp x M x dp x dp x + dp x 




-n 2 


eft 1 dV eft 


dp 2 M 2 dp 2 dp 2 


d 


dp 2 


/■ > 

O 

fin 'hP + ^2 i 
^ 6// 2 y 


-r h 


dp_ 

dp 3 M 3 dp 3 di] 3 dp 3 


1 dV eft d 


fiz 2 r] 2 P + x<h 


dp 


dpj 


= 0 

(4.26) 


With an assumption (Caughey, 1963) 

C fin / MJi ) = ( fin / M 2 <f > 2 ) = (/? 33 / = Y 


equation (4.26) can be solved to obtain the joint probability density of displacements and velocities 
in terms of the transformed coordinate system as 


p(p 1 ,P2,rh,Pifthfth) = ce xP 


-{M x p\ 12 + M 2 p\ I2 + M 3 p 2 3 / 2 + V(p x ,p 2 ,p 3 )\ 
n '■ 


(4-27) 


Perfonning tlie inverse orthonormal transfomiation and noting that the term 


M x p\ / 2 + M 2 tjI I2 + M 3 p] / 2 + V(p x ,p 2 ,p 3 ) 
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on the right hand side of equation (4.27) represents the total kinetic and potential energy of the 
system and that y is a constant, the joint probability density of displacements and velocities in the 

original set of coordinates is 


P(x lt x. 


2 ? 


x 3 ,x x ,x 2 ,x 3 ) 


cexp 


- — {m x x 2 x +m 2 ±l + m 3 x 3 + F(x,,x 2 ,x 3 )j 
n 1 ’ 


(4.28) 


The joint probability functions /?(x,,x 2 ,x 3 ) and p(x x ,x 2 ,x 3 ) are obtained from equation (4.28) 
as (Roberts and Spanos, 1 990) 


P( x i ,x 2 j x 3 ) — j* J J" p(x x ? x 2 , x 3 , x 3 ,x 2 , x 3 ')dx x cfx 2 (ix 3 


= c x exp 

with 

j oo oo oo 

c * = y l!l ex p 


F(Xj , x 2 ,x 3 ) 

71 


V{x l ,x 2 ,x J ) ]fhc l dx 2 dx 3 

TC 


(4.29) 


and 


p(X l ,X 2 ,X 3 ) — j" J" J P( X l-> X 2’ X 3^l’^2^i) C ^ X l^ X 2 ( ^ X 3 


-OO— 00-00 


1 L 

n\2 


^m^n 2 m 3 exp 


'** +W2 ** +m 3 i 3 2 )] 


The joint probability function p(x 15 x 2 ) is obtained as 

CO 

p(x 1>X 2 )= ^ P( x x ,X 2 ,X 3 )dx 3 
—00 

oo r 

= fc,exp {F(Xj,x 2 ,x 3 } dx 3 

J IT 


(4.30) 
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= c 2 exp 


m iY 1 


4 A 


2 _i_^2 _.®13 


— G),, + G) 


■ ^ , '11 1 ^n, 2 | 

n I 2 V G)^,/ 


4 A 


Jl , M ..2 , ^,2 *^23 

Xl + 2 l <y22 -2 


+ 


^..2 

^ 12 2 

V GJ33 ^ 


<y 33 y 

*1*2 +(Gi^A 1 )g(r 1 ) + (G3^A 2 )g(x 2 ) 


with 


c 2 = 


K 


1^33 \ m iY 



exp 


-00—00— 00 


Y_ 

K 




The probability density function p{x l ) and p(x 2 ) are 

00 co 

P(* 1)= j \ pix lt x 2 t x 3 )di 2 di : 3 


’1 



Y f 1 -2] 

-ll 


exp 

— )t w i x i r 

;r v 

2 


. x U J J 


00 00 

Pi* 2 )= j j Pix^x^dx, 


h r 


rfi .2]' 

exp 

| „ WI 2*2 [ 

i 2 

L n l 2 JJ 


The variances of velocity responses x 1 and x 2 are obtained as 

00 

= {*1 P(*i)<&i 


TC 


"hr 

00 

<*\ = jx 2 2 /3(x 2 )rft 2 


tt 


»z 2 r 


(4.31) 


(4.32) 


(4-33) 


(4-34) 


(4.35) 
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Combining equations (4.34) and (4.35) 
n 1 


= 


y 


(4.36) 


Tlie joint probability density function for the displacement responses x } and x 2 , from equations 
(4.3 1) and (4.36), can be written as 


P( x n x 2 ) = c i exp 


1M2 f 1 

a i cr. 2 
■ x \ x 2 V. 


,2 , ^.,2 ® I3 


-- 


V 




x, 2 +- 

2 


V 




CO, 


co^J 


+ 


^ 2 ^ 13 ^ 23 ^ 

*>12 i , 

v <y 33 J 


X V X 1 + (® l * 1 )g(- *1 ) + (® ^2 )g( X 2 ) 


with 


C 2 = 


n 


\& 3 }\ m 3 r - 


oo oo OO 

J I i ex P 


V^2 

L 


■ V(x 1 ,x 2 ,x 3 ) pc ,^ dx 






(4.37) 


4.5 EXTRACTION OF BEARING STIFFNESS PARAMETERS 

The bearing parameters are obtained from the experimentally obtained random response in terms of 
the linear stiffness parameters (Q 1 (Q 2 and the nonlinear stiffness parameters l l , X 2 and the 
disc-bearing mass ratios }J. X and p 2 . These parameters are obtained for both, the vertical and 
horizontal directions, the problem formulation, in the horizontal direction, remaining identical to 
that in die vertical direction. 

The joint probability density function p{x x ,x 2 ) for a set of displacements (x^ ,x 2j ) 

and X, UJ1 >X 1; ). from equation (4.37) are 
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P( x i,,* 2 ,) = e 2 exp 


tJmiMi fl 

L l 2 


ty 


®n +t0 », ~ _ 2 


4 A 
13 


CO 
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x\ ?+- U 2 22 +C» 2 n -Z? 
h 2 v "* a>\ 3 j 


r 2 _(0 2 u a) 2 2 ^ 

^ 12 2 

v < y 


x i, x 2, + «A, )g(x 1( ) + (<y^A 2 )g(x 2 _ ) 


(4.38) 
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(4.39) 


Combining equations (4.39) and (4.38) gives 
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P( x \ ? X 2 ) = 
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(4.40) 

For N values of each displacements, x 1} , x la , ,x^ , and x 2 ,x 2 , ,x 2 ^ , equation 

(4.40) can be expressed as a set of (N - l) 2 linear simultaneous algebraic equations, as. 
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(4-41) 


/ = 1,2, j = 1,2, AN-l) 


Equations (4.41) are solved for Co\ , <o\ it X lt X 2 , and , using the least square fit 

technique. Tire joint probability function, p(x l ,x 2 ) and variances, cr* and a\ are computed 
from the experimentally obtained displacement and velocity data (x l7 x 2 and x 2 ) , which are 
taken as zero mean Gaussian processes and the nonlinear spring force provided by the rolling 
element bearings is taken to be cubic in nature i.e. G(x) — x . The stiffness -matrix defined as for 
the three body lumped-parameter shaft model ( Figure 4. 1) can be defined as (Childs, 1993) 


[K] = m 3 


a>n 


0) 2 n 

CO 2 


31 C °-il lt/ 32 

ro 2 i co 2i co 22 


12 El 


-1 0 

2 -1 

-1 1 


(4.42) 
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4.6 EXPERIMENTATION 

The laboratory rig, described in Chapter 3, is employed for experimental illustration of the 
procedure. Typical displacement and velocity signals, in the vertical direction, picked up (after 
balancing the rotor) by the accelerometers mounted on the two bearing housings are given in 
Figures 4.2, 4.3, 4.4, and 4.5. The joint probability density function, p(x 1 ,x 2 ), of the 
displacements is shown in Figure 4.6. The following set of data is taken for the rotor 

El = 1.03 x 10 8 N-mm 2 

m 2 = 0.515 kg 

L = 250.0 mm.. 

The parameters estimated from equations (4.41) are given in Table 4. 1 . 


Table 4.1 Estimated bearing stiffness and mass parameters 


Parameters 

Vertical 

Horizontal 

a> 2 th (rads/sec) 2 

2.13 x 10 7 

2.12 x 10 7 

A j (mm)' 2 

-1.25 x 10 6 

-1.41 x 10 6 

co (rads/sec) 2 

1.98 x 10 7 

2.21 x 10 7 

X 2 (mm)' 2 

-1.15 x 10 6 

-1.36 x 10 6 


2.10 

2.34 


(Refer to Appendix C for Statistical Error Table) 







Velocity i. (mm/sec X 10') DiSP ‘- *' <m “ X ,0< > 



Figure 4.2 Displacement signal in the vertical direction at bearing 1 



Time t (sec X 10 1 ) 


Figure 4.3 


Velocity signal in the vertical direction at bearing 1 






Time t (sec X 10 J ) 

Figure 4.4 Displacement signal in the vertical direction at bearing 2 



Time t (sec X 10 1 ) 

Velocity signal in the vertical direction at bearing 2 


Figure 4.5 





Joint probab. funct. p{x x ,x 2 ) 



Figure 4.6 


Joint probability density distribution of vertical displacements 
at bearing 1 and bearing 2 
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4.7 MONTE CARLO SIMULATION 


The algorithm is tested by Monte Carlo simulation. The experimentally obtained values of 
if ^2 ■\M 1 M 2 are fed into equation (4.8). Broad band excitation forces, 

/,(Oand f 2 (t) with zero mean and Gaussian probability distribution as described in Figures 
4.7, 4.8, 4.9 and 4. 10, are computationally simulated. These forces are also fed into equation (4.8) 
and the equations are numerically solved through fourth order Runga-Kutta to obtain the 
displacement and velocity responses, x, ,x 2 ,x, and x 2 . These simulated displacement and velocity 
responses, at the two bearings, are shown for the vertical direction in Figures 4. 1 1, 4. 12, 4. 13 and 

4.14. The joint probability distribution of the simulated vertical displacements is shown in Figure 

4.15. The simulated response is now fed into equation (4.41) to obtain the values of 

<y* ,<y^ j,/L 2 and ■Jm 1 ]U 2 . A similar exercise is carried out to obtain the parameters in the 

horizontal direction. These values are listed in Table 4.2. 

Table 4.2 Estimated and simulated bearing stiffness and mass parameters 



Vertical 

Horizontal 

Parameters 

Experimental 

Simulated 

Experimental 

Simulated 

co 2 (rads/sec ) 2 

2.13 x 10 7 

2.12 x 10 7 

2.11 X 10 7 

2.21 X 10 7 

k x (nun )' 2 

-1.25 x 10 6 

-1.41 x 10 6 

-1.29 x 10 6 

-1.45 x 10 s 

oo 2 ^ (rads/sec ) 2 

1.98 x 10 7 

2.21 x 10 7 

1.95 x 10 7 

2.23 x 10 7 

X 2 (nun) ' 2 

-1.15 x 10 s 

-1.36 x 10 s 

-1.63 x 10 6 

-1.86 x 10 6 

1 M 2 

2.10 

2.43 

2.34 

2.50 


The good agreement between the values of the bearing stiffness parameters, ty 2 ,co 2 2 ,A I ,A 2 and 

tJJlju 2 , obtained by processing the experimental data and those from the Monte Carlo 

simulation, indicate the correctness of the experimental and algebraic exercises. It should be noted 
that the simulated values of the bearing stiffness parameters are obtained for an ideal white noise 
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Figure 4.9 Simulated random force at bearing 2 



Figure 4.10 


Probability density distribution of simulated force 
at bearing 2 
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Time t (sec X IQ 1 ) 

Figure 4.13 Simulated displacement signal at bearing 2 



Figure 4.14 


Time t (sec X 10 1 ) 

Simulated velocity signal at bearing 2 




Joint probab. funct. p(x. 



Figure 4.15 Joint probability density distribution 

of simulated vertical displacements 
at bearing 1 and bearing 2 
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excitation, while the experimental ones are obtained by processing the actual response of the 
system, where the unknown excitation was idealized as white noise. It also needs to be pointed out 
that the values of the damping parameters a. ~ , are not required for the estimation procedure 

(equation (4.41) ) . Any convenient set of values of cc ^ , can be employed in equations (4. 8-4.1 1) 
for the purpose of simulation. 

4.8 VALIDATION 

The analytical formulations of Ragulskis, et al. (1974) and Harris (1984), described earlier in 
Chapter 3, are used to validate the stiffness parameters obtained from the experimental signals, 
through the procedure developed. The bearing stiffness expression obtained from the analytical 
formulations given in the literature mentioned above, are listed in Table 4.3, under the column 
"Theoretical Stiffness’. The stiffnesses estimated for the two bearings, through the procedure 
developed, are also listed in Table 4.3. 


Table 4.3 Estimated and theoretical (Ragulskis, et al., 1974; 

Harris, 1984) bearing stiffness parameters 


Preload 

(mm) 

Theoretical Stiffness 

(Radial) 

(N/mm) 

Estimated Stiffness 

(N/mm) 

(at bearing 1) 

Estimated Stiffness 

(N/mm) 

(at bearing 2) 

0.0002 

1.20 x 10 4 - 4.01 x 10 10 x 2 





1.08 x 10 4 - 4.19 x 10'° x 2 

1.01 x 10 4 - 4.91 x 10 10 x 2 



(horizontal) 

(horizontal) 

0.0003 

1.47 x 10 4 - 2.18 x 10 10 x 2 



0.0004 

1.69 x 10 4 - 1.42 x 10 10 x 2 



0.0005 

1.89 x 10 4 - 1.02 x 10 10 x 2 





1.10 x 10 4 - 4.13 x 10 10 x 2 

1.02 x 10 4 - 3.52 x 10 10 x 2 



(vertical) 

(vertical) 

0.0006 

2.08 x 10 4 - 6.09 x 10 9 x 2 
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As emphasized in the earlier chapter, the bearing stiffness is critically dependent on the preloading, 
g, of the balls. The exact value of the preloading of the bearing balls in the shaft-casing assembly, 
especially during operations which have involved wear and tear, would be difficult to determine. 
The stiffness of die test bearing is plotted in Figures 4.16 and 4.17 as a function of the radial 
deformation, x, for various allowable preload values, g. The bearing stiffness obtained 
experimentally, using the procedure developed, also shown in Figures 4.16 and 4.17, shows good 
resemblance to theoretically possible values. It is also to be noted that the theoretical stiffness 
calculations are based on formulations which analyse the bearing in isolation of the shaft. The 
comparison between the experimental and theoretically possible stiffness is also listed in Table 4.3. 
The expressions for die dieoretical stiffness in Table 4.3 have been obtained by curve fitting the 
stiffness values obtained from equation (4.41), through a quadratic in x. 

In the analysis, presented in Chapter 3, where the shaft flexibility has not been accounted for and 
die shaft is treated as a rigid body, the bearing stiffness for the same experimental set-up was 

found to be 1.32 x 10 4 - 5.08 x 10 10 x 2 (N/mm) and 2.23 x 10 4 - 8.50 x 10 10 x 2 (N/mm) in the 
horizontal and vertical directions respectively (Refer Table 3.2). A comparison with the stiffness 
values of Table 4.3 reveals the influence of shaft flexibility. 

While a good agreement on the bearing stiffness parameters is observed between the values 
generated following the method of Harris (1984) and Ragulskis et al. (1974) and those obtained 
experimaitally tiirough die present procedure, die values of the effective masses at the bearing 
ends, obtained as by-products of the present procedure, also look reasonable. The experimentally 

obtained values of are 2.1 and 2.34 (Table 4.1) in the vertical and horizontal directions 

respectively. If die two bearings were taken to be identical, the effective mass computed (knowing 
that die disc mass is 0.515 Kgs) for each of the bearing aids turns out to be 0.245 Kgs in the 
vertical direction and 0.220 Kgs in the horizontal direction. These values look reasonable, in view 
of the fact that along with some contribution from the bearings themselves, a division of the mass 
of the shaft, which in diis case is 0.306 Kgs, is seen at the two bearing ends. 



Stiffness K( Xl ) (N/mm X 10“*) 



Figure 4.16 


Stiffness comparison at bearing 1 

1 2 - Present study (In vertical and horizontal directions respectively) 
3, 4,5, 6, 7 - Harris (1984) and Ragulskis et al. (1974) with prelaod 
0.0002, 0.0003, 0.0004, 0.0005 and 0.0006 nun. respectively (In any 

radial direction) 




Stiffness K(x 2 ) (N/mm X KT*) 



Figure 4.17 Stiffness comparison at bearing 2 

& 1,2 - Present study (In vertical and horizontal directions respectively) 

3*4,5, 6, 7 - Harris (1984) and Ragulskis et al. (1974) with prelaod 
o!(H) 02, 0.0003, 0.0004, 0.0005 and 0.0006 mm. respectively (In any 
radial direction) 
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4.9 REMARKS 

The inverse problem for the flexible rotor, features a nonlinear multi -degree-freedom system and 
the essential feature, of the approach to the inverse problem of parameter estimation is an 
appropriate coordinate transformation, so as to arable the governing equations to be modeled as 
Markov Processes through the Fokker-Planck equations. In addition to the saliart features of the 
estimation procedure already mentioned in the earlier chapter, accounting for the shaft flexibility, 
also enables to obtain estimates of the effective masses at the bearing stations. 



CHAPTER 5 


BEARING STIFFNESS ESTIMATION 
IN MULTI-DISC ROTORS 


The parameter estimation problem in rotors carrying more than one disc on flexible shafts is 
considered next. The governing nonlinear differential equations for such multi-mass flexible 
systems are derived for balanced rotors with random excitation at the bearings. The governing 
equations are subjected to a coordinate transformation and modeled as Markov Processes. General 
form expressions, for the first order probability statistics of the response are obtained. The 
statistical response is processed to extract the rotor-bearing stiffness parameters. 

5.1 EQUATIONS OF MOTION 

A balanced rotor, with n discs mounted on a massless flexible shaft supported in nonlinear 
bearings at ends is shown in Figure 5.1. The shaft is treated as free-free body, carrying unknown 
effective bearing masses and m 2 at its ends and the known disc masses , m 4 , . . .,tn ll+2 ■ 

Following the procedure mentioned in Chapter 4, for incorporating the bearings through external 
“ forces ”, F b , acting on masses m l and m, and taking the shaft parameters to be linear, the 

equations of motion are written as 


[M]{X} +[A){X) + {3' / aX) = {F} 


(5.1) 


where {A} is the displacement vector, [M] and [A] the mass and damping matrices respectively, 

{F} = {F l F 2 0} T (5 ' 2) 

V(x 1 ,x 2 ,---,x (n+2) ) = {H 2 ){X} T {{K\ + [K l \){X} + {5'(^)} r [^.vr3{?( x )} ( 5 - 3 > 

«*)>=■< ([' (J y(W4) m o> T (5.4) 

[K] is the shaft stiffness matrix and the linear and non-linear stiffness matrices of the bearing, 
respectively, are 
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[*J = 

5.2 ORTHONORMAL TRANSFORMATION 

The Markov vector approach extended to nonlinear multi-degree-of-ffeedom systems (Nigam, 
1983) is adopted for the solution of equations (5.1). Equations of motion, (5.1), with damping and 
force F t set to zero, are solved for eigenvalues p 2 , i = l,2,...,(« + 2) and orthonormal modal 
matrix [17] , such that 


0 

o 

o' 

0 

[*vzl = 

k Np 

0 

0 

O' 

0 

0 

0 

' 0_ 


_ 0 

0 

. 0 


[U] T [K\[U] = [p 2 ] 
[U] r [M][U] = [I] 


(5.6) 


where [p 2 ] is the diagonal eigenvalue matrix, while [/] is an identity matrix. 

Application of coordinate transformation 

{X}=[U]{r ] } (5-7) 

and premultiplication by [U] T , the equations of motion, (5.1), yield. 


Vi +PiVi +(l/M)^( 7 7i = r /2 J --%7 ( „ + 2)) / =<h 


i = l,2,...,n + 2 (5.8) 


M, is the modal mass. Modal damping matrix [/?] and the generalized force vector {q} are 
[f3] = [U] T [a][U] 

{q} = [U] T {F} 0 

and the potential energy, in generalized coordinates, can be expressed as 



v (rh,r] 2 ,-, T J(n+ 2 ) ) = (1 / 2){t]} t [U] t ([K] + [K L ])[[/] { 77 } 
+{g([U] { 77 })} r [K m ] {*([17] {;?})} 


(5.10) 
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The approach to obtaining the response of the system is simplified with the engineering assumption 
that the random excitation to the system is such that the generalized forces, q t , in equations (5.8) 

can be treated as ideal white noise. Treating the excitation of equation (5.8) as uncorrelated 
Gaussian, white random forces with the following properties 

£[<7,(0] = ° 

(5.11) 

%, (*1 )?, (*2 )] = 2 ^,£(t 2 - t t ) 


where denotes the excitation intensity factor and S(t 2 -t x ) is a Direc delta function. 


5.3 F-P-K EQUATION 


Rewriting equations (5.8) in state space form as 


drii / dt = hi 

3qja = q, - ft, 77 , - (1 / )(cV / ffrj,) 


i = 1,.2 n + 2 (5.12) 


the drift and diffusion coefficients (refer equations (3.5 - 3.7) ), for equations (5.12) can be written 
as 


a i = hi 

bj = qj - /3 jj hj ~ (1 / Mj )W{q x T ] n+2 ) / drjj ] 

and 

c ,.=0 

d jj =l7t(j)j 

e- =0 
u 


i = 1,.2,...,« + 2 

y' = l,2,...,« + 2 


(5.13) 


z' = 1,.2,...,tj + 2 
j = l,2,...,n + 2 


(5.14) 
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With the help of the above, the joint probability density function, 
p( r h-> T h.-’“' ■> T l( n + 2 ),V\-> :r h,‘",V(n+ 2 )) , f° r the motion governed by equations (5.14) can be 
described by the Fokker-Planck equation 


n+ 2 

z 


1=1 



i gv_ % 

Mi 07] i dT]i 


&ii 


~ \ 


v 


PiiPiP + rf t ~~ 

op i) 



(5.15) 


5.4 RESPONSE 


For a stationary case equation (5.15) reduces to 


n+2 

z 

i=l 


_. dp_ 1 dV_ d 

07] i Mi drti dfji + drj, 


PuhiP + ^i — 

V OT]i 


(5.16) 


With the assumption 

Pu / M\(j) x = P 22 I M 2 (j) 2 = /?(/i+ 2 )(n+ 2 ) ! M {n+ 2) $ (n+2) ~ 7 


equation (5.16) can be solved to obtain the joint probability density of displacements and velocities 
in terms of the transformed coordinate system as 


Pin i? / 72 5 ’"> 7 7 (/i+ 2)5 Hi 


V«)) = cex P ! 


-\z;{p} T [M]{T]}+V{ri x ,T] t 
71 (2 



(5.17) 


Performing the inverse orthonormal transformation and noting that the term 


(l/2){7} r [M]{i7}+F(7 1 ,7 72 ,-,7(„ + 2)) 

in the bracket on the nght hand side of equation (5.17) represents the total energy of the system and 
that y is a constant, the joint probability density of displacements and velocities in the original set 


of coordinates is 



p(x l ,x 2 ,---,x (n+2) ,x l ,x 27 -~,x (:n+2) ) - cexp j-W T [M]{Z}+F(x 1 ,x 2) — ,x (n+2) ) 


(5.18i 


The joint probability functions p(x x ,x 2 ,-,x (n+2) ) and p(x u x 2 ,-,x (n+2) ) are obtained ffc, 
equation (5.18) as (Roberts and Spanos, 1990) 

oo oo oo 

> *2 > ' ' ' > X (n+2) ) = | J”' \ P( X l^ X 2 ^"y X (n+ 2 )^i,X 2 ,-",X (n+2) )dx l di 2 ---dx u+2) 


: i ex P -^{R(Xi,x 2 ,-,.x ( „ +2) )} 


r 1 * J J—Jexp {^(*i,x 2 ,-,x {lI+a) )} cbc x dx 2 -dx 


(5.19) 


p(x x , x 2 , • • • , x ( „ +2) ) = J | | p(x, , x 2 , ■ - , x (n+2) , Xj , x 2 , •• • , x ( „ +2) )dx x dx 2 ■■■cbc l 


exp -U X -{X) T [M){X) 



(5.20) 


The joint probability function p(x,,x 2 ) is obtained as 

oo oo oo 

p(x x ,x 2 ) = J J | p(x x ,x 2 ,---, x (n+2 ) )dx 2 dx 4 • - - dX(, l+ 2 ) 


oo oo oo r y 

1 j J j eX P ~~V( X \’ X 2’'"’ X <,n+2)') dx-sdx 4 •■•^(n+2) 

. -oo — oo — oO - ^ 

: c 2 exp j-(A/ lf22 ! K.- 1122 K + \{ M <n 1 M k\ m) x l + { M xii 1 M k\ m) x i x z 
l rt [2 

+k NL{ g~(x i ) + ^.vr, S ( x 2 )}] 


with 
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It _ 


Vm, 


If 1122 


‘ 00 00 oo 

JM 


exp 


„^( X l’ X 2>'"> X (n+2)) 


JT 


M r 3^4”‘ < ^'(„+2) 


M = [*] + [*d 

and M k denoting the minors of the matrix [at] . 

The probability density function p^) and p(x-,) are 

OO oo OO 

P(X l)= f J-- J p(^i,-^2 ,---,-V (m+ 2))^2^3---^(„+2) 

—CO —CO -oo 



1 

T 

exp 

_r_ < 

fl 41 


1 


T m i x i r 


7t ' 

V 2 


TC 

U JJ 


(5.21) 


(5.22) 


/** 2 )= J f- f/>(JCi;*2,-,i(„ + 2))flSide J -*(» + 2) 

—00 -00 -00 

(5.23) 


1 

It 

exp 

1 

i 

1^ 

>T— ' 
.? 

& 

~ v 1 

K y 

ii 2 


L JJ 


The variances of velocity responses iy and x 2 are obtained as 

00 

= \x\p(x,)dx x 

—OO 

= k / /n,y 

00 * 

C7", = J *2 /?( X 2 )i& ; 

-OO 

= k / /?i 2 y 

Combining equations (5.24) and (5.25) 


= (tt / x )( 1 / V m i" z 2 ) 


(5.24) 


(5.25) 


(5.26) 


the joint probability density function for the displacement responses x,and x 2 , from equations 
(5.21) and (5.26), can be written as 
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P(x l ,x 2 ) = c 2 exp 


<^< 7 * l2 ( Mk22 1 4( jV ^i ' M *iu2 ) x l 

+ (M, I2 / M k1122 )xjX 2 + A vt] g 2 ( Xl ) + A vl _ g 2 (x 2 )}] 


with 




2TCCT k a - 2 Jmpru 


M k 1 122 


=o CO 

j f- Jexp 

— oo-oo — oo 


l ( X l’ X 2>'"’ X (n+2)) 


Vbc z dx,-cbc, n+2) (5.27) 


5.5 EXTRACTION OF BEARING PARAMETERS 


Based on the above analysis, the bearing parameters, namely, the linear stiffiiess parameters 
, k , the nonlinear stiffness parameters Ar iV£j , and the bearing masses and nu , are 

extracted from experimentally obtained random response. These parameters are obtained for both, 
the vertical and horizontal directions. Tire problem formulation, in the horizontal direction, remains 
identical to that in the vertical direction. 


The joint probability density function p(x l ,x 2 ) for a set of displacements (x^ ,x, ) 
(x, ,x,.. '), (x, >x, and x, )x, ) from equation (5.27) are 


p( x^ ,x 2 ) = c 2 exp 


/ — — \'z{Mk 22 I +-{M kU ! 

l 2 2 

+ (M, 12 / M, n22 )x 1( x 2j +k s . Li g 2 (x h ) + k NLz g 2 (x 2j )}] 


(5.28) 


p(x, ,x, ) = c 2 exp 

/ ‘(Ml) ^(,,/tl) 7 * r 


r , AM K n 1 M <1 122 Km, + T( M «'l 1 Mk1 “Km 

<?iV iz y!mx”h >- 2 

+ {J^k\ 2 I M KU2 2) x l u - u X 2 tj . u + k NL] g ( x l„. u ) + ^.VL,S’ ( X 2 (J , tf ) jj 


(5.29) 
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Defining Ax, = x. - x, • Ax, = r _ Y 

for small Ax v and Ax 2 , one can write 


p(x, ,x, ) = c, exp 

r x l d*n 


exp 


M ^ X l + \( M <n/M Kll22 )xl 

+ (M Kl2 / M Kn22 ) Xlt x 2j 
+*VL,r(^)+^V£,.g 2 (^)}] 

{(M k22 /M Kn22 )x h Ax h +(M Kll / M kU 22 )x 2j Ax 2j 


-1 


V ,M 1 OT 2 


+ ( M /cl2 / ^-1122)0^2; + *2,^*1,) 
^.V£.| ^(*1, )^-l, G(X 2 j )Ax 2 ^ } ] 


(5.30) 


Combining equations (5.30) and (5.28) gives 

-1 


P(\. l) ’ x 2 ll .J = P(Xi l >X2,)exp 


z \i.^K22 ^ M kU22 )*[. Ax^ 


(5.31) 


ffi, ^.v, V , "l / ”2 

■*■(^*•11 ^ M Klx 22 )x 2 Ax 2 ^ 

+ (M, n / M v 1122 )(x 1( Ax 2( + x 2 ^ Axj. ) 

+*. VL, G ( X \ >^1, + * Vi» G (*2, )^2, }] 

For A/ - values of each displacement, x 1|5 x u , ,x, , and x 2) ,x 2i , ,x 2f/ , equation (5.31) 

can be expressed as a set of ( N - l) 2 linear simultaneous algebraic equations, as. 



, f p(x h ,x 2 ) 

! — — ln( 

Ax, Ax, p(x, ,x, ) 

h 2 j v-m -un)' 




Kvt , \ 

Axj 
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G(x 2j ) 

Ax, 



{^k22 I M k 1 122 




'x,.> 




iK 


Kll I M Kl 1J2 ) 


C I f *->. 

■ + ■ ' 


_ AXj_. Ax. 






[Ax 2j A x,J 


i = 1 > 1 2 , XN-l) j = 1,2, ,(77-1) 


(5.32) 

Equations (5.32) are solved for , k u , k XL ^ ,k ^ , and -Jrnpnj , using the least square fit 
technique. 


5.6 EXPERIMENTATION 


Two discs are mounted at equidistances on the shaft of the laboratory rotor, as shown in Figures 
5.2 and 5.3. Tire shaft remains supported in tire same bearings that were used for illustration in the 
previous chapters. Tire shaft is driven through a flexible coupling by a motor and the vibration 
signals are picked up (after balancing the rotor) in both, the vertical and horizontal directions, by 
accelerometers mounted on both of die bearing housing. 

The nonlinear spring force provided by die rolling element bearings is taken to be cubic in nature 
(Ragulskis, et al„ 1974) i.e. G(x) = x 3 . The stiffness-matrix of the shaft model of the rig is 
obtained as (Childs, 1993) 

1 0-1 O' 

0 10-1 
-1 0 2-1 

0 - 1-1 2 


[*] = 


\2EI 


El = 1.03 x 10* N-mm 2 


L= 165.0 mm. 


(5.33) 
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Displacement signal in the vertical direction at bearing 1 
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igure 5.6 Displacement signal in the vertical direction at bearing 2 



Time t (sec X 10 1 ) 

Velocity signal in the vertical direction at bearing 2 


Figure 5.7 



Joint probab. funct. p(x 



Figure 5.8 


Joint probability density distribution of vertical displacements 
at bearing 1 and bearing 2 


The disc masses are 

m 3 =0.515 Kg. m 4 = 0.755 Kg. 

(The procedure, however, does not require a knowledge of the disc masses.) 

Typical experimentally obtained displacement and velocity signals (x 1? x 2 ,x,and x 2 ), in the 
vertical direction, picked up by the accelerometer are given in Figures 5.4, 5.5, 5.6, and 5.7. The 
joint probability function, p(x l ,x 2 ) and variances, cr^ and cr^ are computed from the measured 
responses. The joint probability density function, p(x u x 2 ) , of the displacements is shown in 
Figure 5.8. The bearing parameters estimated from equations (5.32) are given in Table 5.1. 

Table 5.1 Estimated bearing stiffness and mass parameters 


Parameters 

Vertical 

Horizontal 

k if (N/mm) 

1.04 x 10 4 

0.87 x 10 4 

. 

. .. 

£. Vii (N/mm 3 ) 

-5.10 x 10 10 

-3.50 x 10 10 

k u (N/mm) 

1.04 x 10 4 

0.86 x 10 4 

k xl, (N/mm 3 ) 

- 3.62 x 10 10 

- 2.19 x 10 10 

(Kg) 

0.21 

0.20 


(Refer to Appendix C for statistical Error Table) 

5.7 MONTE CARLO SIMULATION 

Monte Carlo simulation is employed to check the accuracies involved in making the engineering 
assumptions in the algorithm. The above estimated values of k ^ ,k^ ,k N ^ ,k SL ^ and rnpn^ are 
fed into equation (5.1). Broad band excitation forces, /^(t) and F 2 (?) with zero mean and 
Gaussian probability distribution, as described in Figures 5.9, 5.10, 5.11 and 5.12, are simulated 
on the computer and these are also fed into equation (5.1). Fourth order Runga-Kutta technique is 
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Figure 5.9 Simulated random force at bearing 1 



Figure 5.10 Probability density distribution of simulated force 

at bearing 1 



Probab. funct* p(r 2 ) 


Time t (m-sec) 


.60 


Figure 5.1 1 Simulated random force at bearing 2 
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Figure 5.12 Probability density distribution of simulated force 

at bearing 2 
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Figure 5.13 Simulated displacement signal at bearing 1 
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Figure 5.14 Simulated velocity signal at bearing 1 
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Figure 5.16 









Joint probab. funct. p{x. 



Figure 5.17 Joint probability density distribution 

of simulated vertical displacements 
at bearing 1 and bearing 2 


then employed to numerically compute the displacement and velocity responses, x l ,x 2 ,x l and i. 
The displacement and velocity responses at the two bearings, thus simulated are given in Figures 
5.13, 5.14, 5.15 and 5.16, for the vertical direction. The joint probability distribution of the 
simulated vertical displacements is shown in Figure 5.17. The simulated response is now fed into 
equation (5.32) to obtain the values of k Li ,k I ^,kyi i ,k ;fL ^ and ^m l m 2 . A similar exercise is 
carried out to obtain the parameters in the honzontal direction. These values are listed in Table 5.2. 

Table 5.2 Estimated and simulated bearing stiffness and mass parameters 



Vertical 

Horizontal 

Parameters 

Experimental 

Simulated 

Experimental 

Simulated 

kj (N/mm) 

1.04 x 10 4 

1.49 x 10 4 

0.87 x 10 4 

0.99 x 10 4 

ky^ (N/mm 3 ) 

-5.10 x 10 10 

i 

CO 

CO 

* 

) — * 

o 

-3.50 x 10 10 

-4.81 x 10 l ° 

k u (N/mm) 

1.04 x 10 4 

1.48 x 10 4 

0.86 x 10 4 

0.98 x 10 4 

k yu (N/mm 3 ) 

- 3.62 x 10 10 

-7.83 x 10 10 

- 2.19 x 10 10 

-3.23 x 10 10 

(Kg) 

0.21 

0.26 

0.20 

0.24 


The fairly good agreement between the values of the bearing stiffness parameters, 
k^ , k u , k NL ,k yu and ->Jm x m 2 , obtained by processing the experimental data and those from the 

Monte Carlo simulation, indicate the correctness of the experimental and algebraic exercises. It 
should be noted that the simulated values of the bearing stiffness parameters are obtained for an 
ideal white noise excitation, while the experimental ones are obtained by processing the actual 
response of the system, where the unknown excitation was idealized as white noise. It also needs to 
be pointed out that the values of the damping parameters cty , are not required for the estimation 

procedure ( equation (5.32) ). Any convenient set of values of ccy , can be employed in equations 

(5.1) for the purpose of simulation. 




5.8 VALIDATION 


The bearings in the laboratory rotor rig employed to illustrate the procedure, are the same, as those 
used for illustrations in the previous chapters. Table 5.3 makes a comparison between the stiffness 
parameters obtained by processing the experimental signals from the two-disc rotor with the 
expressions for bearing stiffnesses obtained through the analytical formulations of Ragulskis, et al. 
(1974) and Harris (1984), described in Chapter 3. 


Table 5.3 Estimated and theoretical (Ragulskis, et al., 1974; 

Harris, 1984) bearing stiffness parameters 


Preload 

(mm) 

Theoretical Stiffness 

(Radial) 

(N/mm) 

Estimated Stiffness 

(N/mm) 

(at bearing 1) 

Estimated Stiffness 

(N/mm) 

(at bearing 2) 

0.0002 

1.20 x 10 4 - 4.01 x 10 10 x 2 





0.87 x 10 4 - 3.50 x 10 10 x 2 

0.86 x 10 4 - 2.19 x 10 10 x 2 



(horizontal) 

(horizontal) 

0.0003 

1.47 x 10 4 - 2.18 x 10 10 x 2 



0.0004 

1.69 x 10 4 - 1.42 x 10 10 x 2 



0.0005 

1.89 x 10 4 - 1.02 x 10 10 x 2 





1.04 x 10 4 - 5.10 x 10 10 x 2 

1.04 x 10 4 - 3.62 x 10'° x 2 



(vertical) 

(vertical) 

0.0006 

2.08 x 10 4 - 6.09 x 10 9 x 2 


... 


The theoretical stiffness in the above table is dependent on the preloading, g, of the balls, which 
may be difficult to guess, for bearings which have involved wear and tear under operation. The 
estimated stiffnesses can be seen to lie within the preload range specified by the manufacturer. The 
analytical stiffness of the test bearing is also plotted in Figures 5.18 and 5.19 as a function of the 
radial deformation, x, for various allowable preload values, g. along with the estimated bearing 



Stiffness ^T(x t ) (N/mm X 10 J ) 



Displacement x, (mm X 10 ) 


Figure 5.18 


4, 5,6, 7 - Harris (1984) and Ragulskis et al. (1974) with prelaod 
n 0003. 0.0004, 0.0005 and 0.0006 mm. respectively ( any 


radial direction) 




Stiffness K(x 2 ) (N/mm X 10 4 ) 



Displacement x t (nun X 10 ) 


Figure 5.19 


ffness comparison at bearing 2 . , . 

- Present study (In vertical and horizontal directions respective y) 
' 5,6,7 - Harris (1984) and Ragulskis et aL (1974) with pre laod 
Jui-i n nnnv 0.0004. 0.0005 and 0.0006 nun. respectively (In any 


radial direction) 



stiffnesses. It needs to be emphasized again, that the theoretical stiffiiess calculations are based on 
formulations which analyse the bearing in isolation of the shaft. 

The experimentally obtained values of -J ~rnjrl l , a parameter representing the effective masses at 

the bearing ends, are 0.21 and 0.20 (Table 5.1) in the vertical and horizontal directions 
respectively. If the two bearings were taken to be identical for each of the bearing aids turns out to 
be 0.21 Kgs in the vertical direction and 0.20 Kgs in the horizontal direction. These values look 
reasonable, in view of the fact that along with some contribution from the bearings themselves, a 
division of the mass of the shaft, which in this case is 0.306 Kgs, is seen at the two bearing ends. 

5.9 REMARKS 

Generalized expressions for the governing equations, orthonormal transformation, rotor response 
provide a compact algorithm for parameter extraction. A comparison between the bearing stiffiiess 
values obtained in this chapter, can be made with those obtained in the previous chapters, since the 
test bearings employed in the expenmental rigs are the same. The closeness of results is an 
indication of the correctness of the theoretical and experimental exercises. 



CHAPTER 6 


PARAMETER ESTIMATION IN 

NON-LINEAR ROTOR -BEARING SYSTEMS WITH UNBALANCE 


Estimation of linear and non-linear bearing stiffness parameters, for rotors with unknown unbalance, is 
considered in this chapter. The problem is formulated as single degree of freedom system. The case of a rotor 
with a rigid shaft, in nonlinear flexible bearings is investigated. The excitation to the system consists 
harmonic forces due to the unbalance and random forces due to arbitrary deviations, of bearing contact 
surfaces and subsurfaces, from their ideal design and their progressive deterioration during operation. These 
random forces are comparable to the harmonic excitation forces, if the unbalance is not significantly large. 

The parameter estimation procedure is based on the averaging technique of Bogoliubov and Mitropolsky 
(1961) for deterministic non-linear systems, extended by Stratonovich (1967) for stochastic differential 
equations. The governing equation of motion is transformed from the rapidly varying variables, namely 
displacement and velocity, to variables, amplitude and phase, varying slowly with time. Stochastic averaging 
is done, to take into account the effect of the random excitation multiplied by a correlated term, so as to 
model the slowly varying amplitude as an approximate Markovian process. A second order stochastic 
approximation is carried out and a one-dimensional Fokker-Planck equation is derived to describe the 
Markovian amplitude process. The response to the Fokker-Planck equation is derived and processed for 
parameter estimation. 

6.1 EQUATION OF MOTION 

The analysis is restricted to rotors with the disc mounted on a rigid shaft, the non-linear spring force being 
contributed by the bearings (Figure 6.1). The damping to the system is taken to be linear and the disc carries 
an unbalance Fq . The governing equation for the system is written as 


nix + cx + k L x + k m G{x) = F 0 cos (cot + 0) + y/(t) 


(6T) 




Figure 6.1 


Rotor-bearing model 
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In the equation (6.1), k L and k NL are the linear and non-linear stiffness parameters of the bearings and G(x) 
can be a polynomial in x . The rotor mass is m and co is its rotational speed. The angular location of the 
unbalance with respect to a reference point on the shaft is expressed by 6 . y/(t) is the random excitation 
experienced at the bearing ends. The spring force nonlinearity in rolling element bearings is taken to be cubic 
in nature (Ragulskis, 1974), i.e. 

G(x ) = x 3 (6.2) 

6.2 STANDARD FORM TRANSFORMATION 

The concept of averaging principle, developed by Bogoliubov and Mitropolsky (1961) for deterministic non- 
linear vibration, transforms die equation, involving vibrations which are rapidly varying with time, to a set of 
simple equations for slowly varying response coordinates. This principle, extended by Stratonovich (1967) 
for stochastic differential equations, has been employed to analyse the rotor-bearing system governed by 
equation (6.1). 

Defining 

X = k NL l k L (6.3) 

and since (1/2) is a small quantity (rolling element bearings being highly nonlinear), equation (6.1) can be 
rewritten in terms of the small parameter e = (1 IX) as 

x+co 2 x = sf(x,xX(0) 

where 

/ (x, x, C(t)) = [( F 0 2 / m) cos(cot + 0)+ £(/)2 -2co n X^x-(o) 2 „ -a) 2 )Xx-o) 2 „X x ] ^ 

( 0 2 „=k L tm 
% = c/2ma> n 

C(0 = HO / m 

e being a small quantity, the response can be taken to be harmonic in time with frequency a with slowly 
varying amplitude, A(t) and phase, cp(t) , i.e. 
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x(/) « A{t) cos[<y^ + cp{t)] 

x(t) « -(oA(t)su\[(ot + cp(t)\ (6.5) 

Equation (6.4) can, now, be expressed as a set of standard form equations in terms of the slowly varying 
parameters A(t) and cp(t) as 


A = [x 2 + (x 2 / co 2 )] U2 
cp - - arctan[i / cox] - cot 


( 6 . 6 ) 


or 

A = eG[A,<pX(t)} 
cp = eH[A,<p,£(t)] 


(6.7) 


where 

G[A,cp, £(t)) = G„(A, tp) - {£(t)A 1 0 )}sm{ojt + q>) 

H[A , cp, aOJ = H„(A, (p) - {£(t)Ho)A} cos (cot+<p) (6.8) 

witlr 

G av ( A , cp) = {-((oU / 2) + (a> 2 „ U / 2a ) + (a 2 A 2 A 3 / 4a )} sin 2(cot + cp) 

+ {a „ A/pi } cos 2 (cot +<p)-(F 0 A/ 2ma>) sin(2©t + <p + 0) 

+ (m 2 A 2 A 3 / 8<y) sin 4 (cot + <p)~{(F 0 AI Imco) sin(<p - d) +co „A%A} 

H av ( A, <p) = -(a „A%) sin 2 (cot+cp) + {-(a A I 2 ) +(co 2 A / 2 co) 

+ (a 2 A 2 A 2 / 2r»)} cos 2(cot +cp)-(F 0 AI 2mcaA) cos(2at +<p+0) ^ 

+ (<o 2 ,A 2 A 2 / M cos 4(<y/ + «?) + < 2 ) ~( F o X 1 2mcoA)cos(<p - 0) 

+ (a 2 ,A/2(o)+{ia 2 n A 2 A 2 /Sa)} 
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6.3 STOCHASTIC AVERAGING 

It can be seen that the right hand side of equations (6.7) contain (employing the terminology of (Stratonovrch, 
1967)) ‘oscillatory’ terms, i.e. harmonic functions of cot, along with randomly ‘flu ctuatin g’ terms, i.e. 
— {9(0^ /<y}sin(nrf + <p)and I coA) cos(cot + <p) , which contain the random force term £(t) . 

However, due to tire presence of e in the equation (6.4), the parameters A and cp vary slowly with time and 
can be assumed to remain constant over a cycle of oscillation. The averaging process for A and cp can be 
carried out in two stages - by stochastic averaging and elimination of the randomly fluctuating terms 
involving random force term C(t ) and then averaging over a cycle of oscillation for removal the oscillating 
terms involving harmonic functions of cot . 

The approach to obtaining the response of the system, including stochastic averaging, can be simplified by 
providing arguniarts similar to those in the previous chapters and treating the random excitation to the 
system as ideal white noise with zero mean and Gaussian distribution. 

For a zero mean random excitation C(t ) , the expressions (6.7) for A and cp can be stochastically averaged 
to write the non- ’fluctuating’ amplitude term, A nJ and the non-’fluctuating’ phase term , <p nf as 


<P nf = eH m (A,<p) 

Further carrying out the averaging process over a cycle one gets, the non-’fluctuating - non- oscillatory 
amplitude and phase terms as 

2rt In) 

Kf-no 1 

4 0 

SV-K F 0 / 2mco) si n(p - 0) + (co „£4)j 

2rt In) 

IK 0 

= [-(m / 2)-(F 0 / 2nnoA) cos(<p - 0) +(o) 2 „ /2o))+Oco 2 n X A l&o)] 


( 6 . 11 ) 
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Considering eqnattons (6.7), (6.8) and (6.9) and putting . fluctuatlng , temls back mt<> ^ above 
express, ons, the nen-'oscillating’ approximations for the amplitude and phase are 

A„ o =-[(F 0 ^in(o)M<p-6)Hco,M)]-{at)/oj}sin(a}t+<p) 

V,w=[-(«>/ 2 )-{F 0 / 2 mcoA)cos(<p- 0 ) + {col/2co) + (2co 2 n A A 1 IZco)] ( 6 . 12 ) 

I coA}cx>$, (a>t+<p) 

The equations (6.12) can be condensed as 

A IW = ^\[A,(p\- {^t) I o))sm{G)t +<p) 

<P„o = \. a AP\ ~ {C(0 / (oA) cos(cot + (p) 

witli (6.13) 

G, [A, <P\ = -l(^'„A/ 2 ma>) sin(p - 0) + a> „A£4] 

H x [A,(p J = / 2) - (F„ A / 2nuoA) cos {<p -0)+(co 2 n A/ 2 co) +(3a> 2 A 7 A 2 / 8<y)] 

Tlie ‘truncated’ equations (6.11), giving A nf _ no and (p nf _ no or the ‘truncated’ equations (6.12), giving 

A m and <p n0 can be taken as approximations of A and cp . However, either approximation does not, 
adequately, reflect the influence of non- linearity in the system, for while the expression for (or (p m ) 

does involve the non-linearity parameter A, the one for A n j_ no (orA no ) does not. A higher order of 
approximation for A and ip is, therefore, essential to adequately represent the effects of the non-linearity on 
the statistical characteristics of the response. 

6.4 SECOND ORDER AVERAGING 

For a higher order approximation, instead of obtaining A no and ip no as in equation (6.12), the procedure of 
asymptotic method developed by Bogoliubov (Stratonovich, 1967) can be employed, whereby the non- 
’fluctuating" approximations A nf and ip nf , of equation (6.10) are taken to have the form 
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A nf =A*+eu{A\<p*) 

(p„f =<P* +ev(A\<p*) 

where A * and <p * are expressed as 

A* = 6 ^l(A\<P*) = 6iJ*(A*,<p*) + £ 2 G*(A\(p t )+ ... 

V ~ &H m{A , <p ) = sH\ {A , (p ) + s 1 H 2 (A ,(p* ) + 

Similarly die variations u and v are expressed in series form as 
u{A\(p*) = u x {A\<p') + *h{A\<p*)+ ... 
v( A \(p*) = v, (A \(p*) + tv 2 (A* ,<p*)+ ... 


(6.14) 


(6.15) 


(6.16) 


Equations (6.14), along witli the series expansions of equations (6.15) and (6.16) are substituted into the 
stochastically averaged equations (6.10). 


Noting that 

u = (<•!</ (H* )A* +(rl/ / <\p* )((o + (p * ) 


v = (cK> I (M * )/i* +{&> I (T(p* )((o +<p*) 


(6.17) 


the terms with identical power of e are equated to obtain the equations governing the successive 
approximations. The equation governing the terms involving the first approximation (of order e l ) is, thus, 
obtained as 


G\(A* ,<p*)+(t)(cht x / ckp*) = G* m {A* ,(p*) 

( 6.18) 

H\{A* ,<p* ) + a>{cb , / dp*) = H* av (A t ,(p*) 

(In tire above equation, tlie term G m ,(A,<p ) , of equation (6.10) has been transformed to G av (A ,<p ) ) 

The right hand side of equations (6.18) involve G m and H m , which contain both, the oscillatory and the 
‘non-oscillatory’ terms. The functions u{A\<p*) and v{A*,(p*) are now chosen in such a way that 
G;<,A\<p’) and H*(A* ,<p*) contain no oscillatory terms. Thus, terms involving {&/#?) and 
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(fr/Op*) are equated to the oscillatory parts of G* m and H*„ , and G* X {A* ,cp*) and H[ 

equated to the non-oscillatory parts, to obtain 

G* x (A*,<p*) = -[(F 0 A/2/nG>)sin(p* -6)+co n A$ 4*] 

Hi {A * , tp* ) = 1 / 2) - (F 0 A / ImcoA * ) cos(p* - 0) +(<y 2 A / 2<y) + (3cy 2 A 2 >4 * 2 / 8ry)| 

and 

fy($q(/l ,<p ) / dtp ) = {—(a>AA /2) + (<y 2 A/4 /2ry)+(<y 2 A 2 j4* /4ft>)}sin2(ry/ + p*) 

+{<y „ A£4 * } cos 2(<yt + ) - (F 0 A / 2/nry) sin(2r»t + tp* + 9) 

+(cuf 1 A 2 A 3 /8co)sin4(tvt + <p) 

co{ck) x (A* ,tp*) I dtp* ) = ~(co „ A£) sin 2(<yf + <p* ) + {-(<yA / 2) + (co 2 A / 2<y) 

+(<y 2 A 2 /l* 2 / 2<y)} cos 2(ryt + p * ) - (F 0 A / InuoA * ) cos(2&>/ + + 0) 

+(<y 2 A 2 /I * 2 / 8&>) cos 4(ryf + </>*) 

Comparison of equations (6.19) and (6.13) reveals 

G* l (A*,tp*) = G l (A,tp) 

H* x (A*,tp*) = H l {A,tp) 

Equations (6.20) give 

u x (A*,tp*) = {{M* / 4) -(colM* / 4<y 2 )-(<y*A 2 A* 3 / 8co 2 )} cos2(cot + tp* ) 

+(<y „ A£4 *7 2co) sin 2(cot +<p*)+(F 0 A/2nuo 2 ) cos(2 cot + <p* +9) 

-{co]A 2 A* 3 / 32<y 2 )cos4(<yt + <p*) 


\<p*) are 


( 6 . 19 ) 


( 6 . 20 ) 


( 6 . 21 ) 


( 6 . 22 ) 
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V! (A * , <p * ) = {o)„X% / 2co) cos +p * ) + {-(A / 4) +(® * A. / 4r» 2 ) 

+(<y n A ^ /4<y 2 )}sin2(<yr+i}o )-(F 0 A/2ma> 2 A*) sin(2<y/ + <p* + 0) 

+(o) 2 A 2 A / 32<z> 2 )sin4(cot + tp*) 

For second approximation (of order e 2 ), the governing equations are, similarly, obtained as 

G* +ft>(<3/ 2 tty*) = (#C / dA*)u y +(dG* m /<3p>, -(<3i, / <34*)G* -{duj d<p')H\ 

H *2 +<y(^V 2 I ckp ) — (cH av / <34 )u, +{cH m / <%> )v, -($>, / r34*)G* -(c^, I d(p*)H* 

from which the non-oscillatory G’(A‘,<p‘) and //* ( A \ * ) are obtained as 

Gj{A* ,(p*) = (F' o A/4/ ? ny 2 )H5cyA/4) + (5cy 2 A/4©) + (i* 2 a>^ 2 / < y)}sin(^*-0) 
+(5ry (1 A£F 0 / 8 mco 2 ) cos(<p* - 0) 

H* 2 (A\<p*)= {(3A*V 2 A 3 /%a))-(3A t2 coy /8<y 3 )-(5l/ 4 ^ 4 / 256ro 3 ) - (ryA 2 / 8) 
-(ft;; 1 , A 2 /8ffl 3 )}-{(F 0 A 2 /4/nru4*)-(F 0 <y 2 A 2 /4wry 3 A*) 

-(17A*F 0 <y 2 A 3 / 32wru 3 )} cos(</>* - (?) - (F 0 ru„A 2 £/ 2mco 2 A * ) sin(^* -0) 

Substitution from equations (6.19) and (6.25) into equation (6.15) gives 

A* = -4(F 0 A / 2 wry) sin(^’ - 0) + ry „A£4* ] + £ 2 [(F 0 A / 4wro 2 ){-(5©A / 4) 

+(5<y 2 A/4ry)+(A* 2 <y 2 A 2 / <y)} sin(p* - 0) +(5ry, ,A 2 £F 0 /8mry 2 )cos(<p -0)] 

= 6'[-(ryA/2)-(F 0 A/2wry4*)cos(<p* -0)+(t» 2 A/2o))+(3ry 2 A 2 .4 /8ry)] 

+£ 2 [{(3A* 2 ru 2 A 3 /8ru)-(3A* 2 ru 4 A 3 /8ru 3 )-(51A* 4 ru^A 4 /256<y 3 )-(ruA 2 /8) 
-(ro^A 2 /8ry 3 )}-((F 0 A 2 /4w<yA*)-(F 0 ry 2 A 2 /4mryV) 

-(17A*F 0 ry 2 A 3 / 32wry 3 )} cos(/ -0)-(F o <y„A 2 £/2w<fl 2 A )sm(^ -0)] 


(6.23) 


(6.24) 


(6.25) 


(6.26) 


(6.27) 
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Noting equation (6.14) and that the ‘oscillatory’ terms are confined to the variables u and v, the non- 
’fluctuating ‘ - non-’oscillatory’ amplitude term A nf _ no becomes 


where A * is given by equation (6.26). 


Consideration of equations (6.28), (6.7) and (6.8) enables writing the non-’oscillatory’ amplitude term, A no , 
as 

A no = Anr-«o -e[(C(t)Alco)sm(cot+<p*)\ 

= -4(^0 AI2nuo) sin(<p* -0)+a> n A^A *] + £ 2 [(F Q JU Amco 2 ){-(5coA. / 4) (6 2g 

+(5co 2 /L/ Aco)+(A co 2 A 2 / <y)}sin(ip* -0)+(5co„A 2 ^F o / 8ma> 2 ) cos(<p* -0)] 

- 4(4VM I (*>) sin(<yf + (p*)] 


Similarly 

(p„o = s[< 0 )A / 2) ~(F 0 A / 2/»r»4*)cos(p* - 6 ) + (co 2 n A/2co)+(3co 2 n A 2 A* 2 / 8 co)] 

+£ 2 [{(3 A* 2 (O 2 , A 2 /8co)-(3A* 2 co 4 X 3 l8co 3 )-(5\A* 4 co 4 n X 4 1 256co 3 )-((oA 2 / 8 ) 

-(co 4 A 2 /8rr) 3 )}-{(F 0 A 2 / AmcoA*)-(F Q co 2 „X 2 I Arm 1 A*) (630) 

~(\1 A* F 0 (O 2 n A 3 / 32/n<y 3 )} cos(<p* -6)-(F Q a n A 2 $l2mco 2 A')$m((p* -6)) 

-e[{^(t)A t coA*} cos (cot + cp * )] 

The expressions in equations (6.29) and (6.30) are taken as approximations of the amplitude and phase 

terms, A and cp , i.e. 

A * A "" (6.31) 


<P *<P,,o 
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6.5 F-P-K equation 

The amplitude term, A and phase term, cp approximations in equation (6.31) are correlated with random 
excitation force C,(t ) . However , since £(t) is assumed to be broad band random process, its correlation time 
is much smaller than the time constant characterising the rate of change of amplitude A* and phase <p , which 
are slowly varying functions of time. It can be assumed that the values of 0) are statistically independent of 
the values of A* i.e. amplitude A' can be approximated as a Markov process (Stratonovich, 1967; 
Khasminiskii, 1966; Papanicolaou and Kohler, 1974 and Lin, 1986). In addition amplitude changes much 
more rapidly than the phase, and hence the amplitude manages to establish an equilibrium amplitude 

distribution p(A*\(p*) for every value of phase <p. 


The modified FPK equation for the simplified equation ( of the form as equations (6.29-6.31) ) of the 
approximate Markov process is given as 


ep(A\tp')/<* = -01 04 * {(«, )p) - 0/ dp *{(£,)} + £’ 2 d 2 / OA * 2 {(c, , )p) 
+t: 2 0 2 / OA* Op* {(d ,,) p) + e 2 0 2 / ckp* {(*„)/>} 


(6.32) 


where modified drift coefficients are expressed as 


a =[sG;+<t 2 G 2 *+...] + £ 2 \K[m* im\p t }dr + £ 2 ] K[OA* / Op* ,^* Z W 
b x = [eH* + e 2 H*+...] + e 2 \K[Op* IdA* ,A*]dt + s 2 \K[d(p I Op ,<p z W 


(6.33) 


and modified diffusion coefficients are 
c u =t' 2 J A'l/i * , A * ]r/r 

-oo 

d u = s 2 | {K[A* ,(p* z ] + K[<p ,A z ]}dt 

— oO 

e u = s 2 \K[(p\(p\]dx 


(6.34) 


—OO 



99 


In equations (6.33) and (6.34) r ( = fe - *,) ) is the time shift and K[A\<p' T ) is the correlation ftmction. 
While computing the correlation function, tire randomly fluctuating parts of the expressions A* and <p' 
(equations (6.29-6.31) ) , i.e. the terms ~{^)XI (o}%m{cot + (p ) and I cA} cos{aX + <p) , which 

contain the random force term 4X0 , need to be accounted for. 

For the one dimensional approximate Markov process , A* of equations (6.29) and (6.31), the drift 
coefficients are obtained as 

a x = {-e(F a XI liiuo) sin(<p* -())- e( ( o„X%A*) + e 2 (F a XI Anno 2 ){-{5coX I 4) 

+{5(o 2 n AI Aw) +{ A* 2 (o 2 n X 2 I a)} sin (<p* - 6) + s 2 {5co n A 2 £F Q / 8 mco 2 ) cos(<p* - 6) (6.35) 

+i- 2 {S(£;co)X 2 !%(o 2 A')) 

b x =0 

The diffusion coefficients are computed to be 

Cjj =e 2 {S(C(0)A 2 /8<w 2 } 

d n = 0 (6.36) 

e n = ® 

where the spectral density of noise 4X0 at the frequency co , defined as 

00 

S(at),co) = 2 J< auKih) > cos<y(t 2 -t y )d(t 2 -f,) (6.37) 


Tire Fokker Planck equation for amplitude A\ from above equations , can now be written as 


-(< 9! dA * )[{-£( F„ A / 2 mco) sin (<p* -0)- s{co n X^A * ) + e 2 (F 0 X / 4 mco 2 ){-(5 coX / 4) 
+(5<y 2 XI A«>) + ( A* 2 (0 2 „X 2 I co)} sin (/ -0) + s 2 (5co „X 2 £F 0 / W ) cos(/ - 6 ) 
+e 2 {S{£ <o)X 2 / %( 0 2 A* }} p] + e 2 {S(£ (o)X 2 I %co 2 ){d 2 p 1 3A* 2 ) = cp(A*)l fit 


(6.38) 
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6.6 RESPONSE 

For a stationary case equation (6.38) reduces to 


-(<?/ dA )[{-^P' a A,l2nuo)^m( < (p* -0)-£^o n ^*)+s 2 {F 0 XIAmm 2 ){-{5coXIA) 

+(5m \X / 4m) + ( A * 2 m 2 n X 2 / m )} sin(<^ * - 0) + s 2 (5m n £ %F 0 / 8mm : 2 ) cos(p * - 0) 

+ £ 2 {S(C; co)X 2 /8m 2 A*}} p] + s 2 {S(^; co)X 2 / 8m 2 }(<9 2 p / M * 2 ) = 0 

The solution to the stationary Fokker-Planck equation, (6.39), turns out as 

p(A*) = cA exp[-{8m 2 / e 2 X 2 S(£,co)}{s{A* F 0 A/2mm)sin(<p* -0) 

+e{A* 2 a) „X£, / 2) - e 2 (FqX I Amco 2 ){-(5^4 * coX / 4) + (5 A *co 2 X / 4m) 

+(A* i co 2 ,X 2 / 3m)} sin(^>* -0)-e 2 (5A*q> n X 2 %F 0 /8wm 2 )cos(ip* -#)}] 

6.7 PARAMETER ESTIMATION 

Tire probability density function for any two values A* and A* +1 ,of the amplitude (with A* +J 
written from equation (6.40), as 

p(A*) = cA* exp l~ {8m 2 / £ 2 A 2 S(£,co)}{tiA{F 0 X/2mco)sin(<p* -0) 

+£(A* 2 m„X^/ 2)-£ 2 (F 0 X/ Amco 2 ){—(5A i coX/ 4) + (5/1,. co 2 „XI Aco) 

+(A* 3 co 2 X 2 /3©)}sinfo>* -0)-e 2 (5A*co„X 2 4F o /8mm 2 )cos(<p* -0)}] 

p(A* +l ) = cA* +l exp [-{8m 2 / e 2 X 2 S(£,co)}{e(Al l F 0 X/2mco)M<P* -0) 

+s( A ; 2 co „ X%/2)- e 2 (F 0 X I Amco 2 ) {-(5/C, coX / 4) + (5A* +1 co 2 „X / Aco) 
+{AZco 2 n X 2 / 3m)} sin(^* - 0) - e 2 (5A- +l co„X 2 ^F 0 / 8mco 2 ) cos(/ -0)}} 


(6.39) 


(6.40) 


) A*) , can be 


(6.41) 


(6.42) 
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Defining, AA, ( A M A ,. ) , for small A A* , one can write, from equations (6.41) and (6.42) 

[M< 1 )/M<)J = (^,^;)exp[-{8ru 2 / £ 2 A 2 1 S(4:,ru)}{ f (A4;F 0 A/2 WW )s.n(^*-t?) 

+s( A, AA, co „ . A£) - e 2 ( F 0 2 / 4/n<y 2 ) {-(5 A4 * col/ 4) 

, j (6.43) 

+(5A^, (0„A/4(o)+(A* AA](o\£ / <y)} sin(<p* - 0) 

-s 2 (5AA i <y„2, 2 ^F 0 1 8mco 2 ) cos{cp* -<?)}] 

For N amplitude values, A x ,A 2 ,..., A N , equation (6.43) is expressed as set of (A^ — 1) linear 
simultaneous algebraic equations, 

[( 1 / 8A A* co 2 ) In { A; +l p( A*) I A ; p( )} ] {S(C, m) / co „ %F 0 cos 0} 

+[ 1 3 / 1 6w<y ] { sin(<jo * - 0) / <y „ £ cos <?} - [5 / 1 6»i6r 3 ] {co n sin {(p * - 9) / £ cos 0} 

(6.44) 

-[A* 2 /Anno 2 j{n>„Asin(^* -0 )/£cos0}+[aI*]{1/F o costf} +[5sin<p* /8m<y 2 ]{tan6>} 

= -[5 cos <p* 1 8m o ) 2 ] i = 1,2, . . . ,(A - 1) 

Equations (6.44) are to estimate the parameters co n , A, F 0 , £, S(£,(o) and 6., using the Least 
Squares procedure. 

The amplitude displacement and velocity data (x and x) is obtained experimentally and using equation 
(6.6), the amplitude A , and phase (p are computed along with the probability function, p(A ) , to be fed into 

equation (6.44) for parameter estimation. However, equation (6.44) involves A*, (p* and p(A') and as an 
initial approximation the experimentally obtained A, <p and p(A) are taken as A* , <p* and p(A") 
respectively. 

6.8 EXPERIMENTATION 


The procedure described above is illustrated on the laboratory rotor. A disc of mass m— 0.41 Kgs. is 
centrally mounted on the shaft. The shaft, as in the previous illustrations, is mounted in the same bearings at 
its ends. Accelerometers are mounted on one of the bearing caps to sense the vibration signals. A reference 
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signal is picked up from the shaft by a non-contact eHHv ■ . , _ 

& y on comact Way current proximity probe. The rotor is dynamically 

balanced and then a known unbalance mass is atiached at the disc. The rotor configuration is same as in 

Chapter 3. -The shaft is rotated at a particular speed and the signals are picked up. The eapenmen, is 

repeated for different set of known unbalance masses and for different speeds. 

Experimentally obtained displacement and velocity sigmtls along with the corresponding reference signal are 
shown, for a rotor speed of 1800 rpm, in Figure 6.2. 

6.8.1 Initial Approximation 

The amplitude A and phase cp signals computed from these measured data are shown in Figures 6.3 and 6.4. 
The probability density function, p(A ), of the amplitude is shown in Figure 6.5. The bearing parameters 
estimated, for initial approximation of using experimentally obtained A(t) asA*(t) and cp(t) as (p\t) in 
equations (6.44) are given in Table 6. 1 . 

Table 6.1 Experimentally estimated parameters (after initial approximation) 




Parameters estimate! 


■ 

in e 

gin~citi 

(at 0 = 0°) 

<P 

degrees 

0 

degrees 

m e 

gm-cm 

$ 

k(x) 

N/mm 

1800 

10.5 

0.7 

0.1 

16.80 

0.062 

1.21 x 10 4 - 0.58 x 10 lo x 2 

- 

17.5 

1.4 

0.2 

22.73 

0.071 

1.23 x 10 4 - 0.83 x 10 lo x 2 

- 

24.5 

6.9 

0.1 

30.64 

0.022 

1.39 x 10 4 - 1.06 x 10 10 x 2 

1400 

10.5 

3.0 

0.3 

17.24 

0.083 

0.86 xlO 4 - 0.84 xl0 10 x 2 

- 

17.5 

3.3 

0.5 

23.45 

0.082 



24.5 

4.6 

0.2 

16.56 

0.042 

1.12 x 10 4 - 0.63 x 10 10 x 2 










































Phase (p (degrees X 10 1 ) 


Figure 6.3 


Time t (sec X 10 1 ) 

Amplitude variation of the measured response 



Figure 6.4 Phase variation of the measured response 





Probab. funct. p(A) 
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6.8.2 Iterative Scheme 

For a more accurate estimation of estimated parameters an iterative scheme is devised. The parameters 
estimated from the above initial approximation are substituted in equations (6.22, 6.23), along with the initial 
assumption of taking the experimental A and (p to be A and cp * respectively, to compute the 
variations u and v. These values of u and v and the experimental A and cp are employed in equation 
(6.14) to get a new approximations for A* and cp' ■ The new approximations. A* and <p* , are now 
employed in equation (6.44) for a fresh parameter estimation and thus, the iterative cycle can be continued. 
The final set of parameters estimated, after such iteration is given in Table 6.2. 

Table 6.2 Experimentally estimated parameters (after iteration) 


Speed 

Unbalance 

Parame 

ters estim 

ated 

CO 

rpnt 

m e 

gm-em 

(at 0 = 0°) 

<P 

degrees 

0 

degrees 

m e 

gm-cm 

§ 

k(x) 

N/mm 

1800 

10.5 

2.4 

0.3 

14.23 

0.023 

1.31 xlO 4 - 0.92 x 10 10 x 2 

- 

17.5 

1.6 


20.32 

0.036 

1.41 x 10 4 - 1.10 x 10 10 x 2 

- 

24.5 

3.2 

0.2 

27.32 

0.032 

1.60 x 10 4 - 1.28 x 10'°x : 

1400 

10.5 

1.7 

0.2 

15.43 

0.063 

0.94 x 10 4 - 0.94 x 10 10 x 2 

- 

17.5 

3.2 

0.4 

21.40 

0.063 

0.93 x 10 4 - 1.01 x 10 10 x 2 

- 

24.5 

1.8 

0.1 

23.56 

0.037 

1.37 x 10 4 - 0.72 x 10 10 x 2 


'( Refer to Appendix C for Statistical Error Table) 


6.9 VALIDATION 

The closeness, between the known unbalance introduced in the totor and its angular location and tire 
corresponding oxporimentallv estimated values, is a measure of the validity of the procedure de op 
estimated damping ratios also appear physically reasonable. The beanng stiffness parameters obtamed in tire 
prevrous Chapter 3, employing tire analytical formulations of Ragulslds et al. (1974) and Hams (1984), are 

reproduced in Table 6.3, for easy reference. 
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1 able 6.3 Theoretical (Ragulskis et al., 1974 and Harris, 1984) 
bearing stiffness parameters 


Preload 

(mm) 

Theoretical Stiffness 

(Radial) 

(N/mm) 

0.0002 

1.20 x 10 4 - 4.01 x 10 10 x 2 

0.0003 

1.47 x 10 4 - 2.18 x 10 10 x 2 

0.0004 

1.69 x 10 4 - 1.42 x 10 10 x 2 

0.0005 

1.89 x 10 4 - 1.02 x 10 10 x 2 1 

0.0006 

2.08 x 10- 6.09 x 10 9 x 2 


The estimated bearing stiffnesses of Table 6.3 can be seen to show good resemblance to theoretically 
possible values of Table 6.3. 

6.10 REMARKS 

The analysis presented underlines that stochastic averaging techniques can be usefully employed to address 
the inverse problem of parameter estimation in a rotor system governed by a nonlinear equation having both, 
harmonic and random forms of excitation. Also, along with the bearing stiffness parameters, estimates of the 
unknown unbalance of the rotor, its angular location and the damping ratio are also obtained, as by-products 
of die procedure developed. 




CHAPTER 7 


CONCLUSION 

Tbe P reseilt stl,d y has concerned itself with the inverse problem of parameter estimation in nonlinear rotor 
bearing systems from the vrbration signals measurable at the bearing caps. Methods in statistical dynamics 
have provided tire framework of analysis. Nonlinear rotor bearing system have been modeled as Markovian 
processes or as approximate Markovian processes after stochastic averaging to render standard form 
equations, which were solved to obtain the response statistics. 

Tire following rotor systems were considered 

(i) Balanced Rigid Rotor 

(ii) Balanced Rigid Rotor With Flexible Disc 

(iii) Balanced Multi-Disc Flexible Rotors 

(iv) Unbalanced Rigid Rotor. 

The parameter estimation procedures, were illustrated on a laboratory rotor model and the results were 
validated through comparison with those from available analytical guidelines. A Monte-Carlo simulation was 
also done to check tire accuracies of the assumptions and the approximations involved. 

In the first three cases, balanced rotors (or rotors with negligible unbalance) were considered and estimates of 
the nonlinear stiffness parameters of the bearings were obtained along with the estimates of linear 
parameters. Tire appropriateness of statistical methods and good estimates of the parameters, provided the 
motivation for taking up the case of a rotor with an unknown small, but not negligible unbalance (If the 
unbalance is noticeably high, then tire priority is to balance the rotor). The statistical averaging technique and 
subsequent iterative refinements, adopted for this case, yield not only good estimates for the linear and 
nonlinear stiffness parameters of tire bearings, but simultaneously also provide good approximations of the 
magnitude of tire unknown unbalance as well as its location with respect to a reference point on the shaft. 
Reasonable values are obtained for the damping ratio too, as by-products. However these damping ratios 
need to be validated through an independent set. 
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These aspects highlight the suitability of the procedure, for usage as an on-line tool for parameter estimation, 
in contrast to other available techniques which either need the machine to be stopped to give it a known 
excitation or others which provide analytically possible ranges of stiffness of the bearing treated in isolation 
of the shaft. The possibility of monitoring the nonlinear part of the bearing stiffness for incipient bearing 
faults and its wear and tear needs to be explored further. 

The unbalanced rotor analysis has been restricted to the case of a rigid rotor in flexible bearings, which could 
be treated as a single-degree-of-freedom system. The analysis needs to be extended to multi-degree-of 
freedom systems, by considering rotors on flexible shafts carrying more than one disc. If such an extension of 
the present scheme can be made, then along with the stiffness estimates, the possibility of obtaining 
approximations for the imbalances on the various discs (as obtained in the present study for a rigid rotor) 
need to be investigated. 

Further, the problems have been formulated by taking the equivalent form linear viscous damping. The role 
of other forms damping in tire analysis needs exploration. The present analysis, based on processing random 
vibrations, is restricted to rotors supported in rolling element bearings. The suitability of the analysis, for 
rotors in fluid film bearings needs to be investigated. 

Instability phenomenon in rotors experiencing random excitations, and the influences of bearing stiffness 
nonlinearity have not been investigated in the present study. However, the present study, and its possible 
future extensions, can facilitate such investigations, by providing estimates of the parameters which form the 
causes of rotor instabilities. 
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appendix a 


INSTRUMENTATION 

A.1 SPECIFICATIONS OF INSTRUMENT: 


I. Accelerometer (Piezo-electric): 

Make - Bruel & Kjaer, Denmark 
Charge Sensitivity - 10.0 pC/ms" 2 
Construction - Delta Shear 


3. Digital 'radiometer: 

Make - Teclock Corporation, Japan 

4. Phase Meter: 

Make - Bruel & Kjaer, Denmark 
Inputs - Two BNC Sockets - 

Phase Accuracy - 

5. Tunable Baud Pass Filter: 

Make - Bruel & Kjaer, Denmark 
Input - BNC Socket 
Selective Bandwidth - 


Type - 4370 

Frequency Range - 0. 1 Hz. ~ 5.0 kHz. 
Maximum Acceleration - 0.6 km/sec 2 

India Type - Eddy current 
2 mm peak to peak. 

50 Hz to 1 kHz 

Range - 0 to 10 000 ipm 

Type -2977 
A - Reference input 
B - Unknown Phase input 
+0. 1° absolute error. 

Type - 1621 

Frequency Range - 0.2 Hz. ~ 20 kHz. 
3% or 23% (1/3 Octave) 


2. Proximity Probe: 

Make - National Aeronautical Laboratory, 
Maximum Amplitude Allowable - 
Frequency Range - 
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6. Charge Amplifier: 

Make - Brucl & Kjaer, Denmark 
Calibrated Output ratings - 


Frequency Range - 


7. Photoelectric tachometer Probe: 

Make - Bruel & Kjaer, Denmark 
Transducer Type - 

Minimum Sensitive - 
Response Time - 

Connector - 

8. Digital Storage Oscilloscope (DSO): 

Make -Larsen & Toubro and Gould Electronics 
Input - 2 Channel BNC Socket 
Screen Update Rate - 30 Traces per sec. 
Interface Bus - IEEE-488 

9. GPIIJ-PC: 

Type - IEEE-488 

Maximum Frequency Range - 100 MHz. 
Transfer Rate - 


Type - 2635 

Selectable in 10 dB steps, 

0.1 mV to 1 V/ms' 2 (Acceleration), 

10 mV to 100 V/ms" 1 (Velocity), 

0.1 mV to 10 V/mm (Displacement) 

Switchable 

0.2 or 2 Hz. to 100 kHz.(Acceleration), 

1 or 10 Hz to 10 kHz (Velocity) 

1 or 10 Hz to 1 kHz (Displacement) 

Type - MM 0012 

Combined infra-red light source and photo-sensor, 
fitted with infra-red filter. 

100 mV at 10 mm distance from a flat white card. 
200 ps minimum for full output. Equivalent to a 
10 mm long reflecting surface passing at 50 m/sec 
BNC with double shield (BNT) 

Type - 4072 

Frequency Range - 0 ~ 100 MHz (DC) 

Maxi. Sampling Rate - 400 M samples/sec 


Connector - Shielded 24 pin conductor cable 
Connected to - PC/AT 
5K bytes/sec in binary mode 



APPENDIX B 


EIGEN VALUES AND EIGEN VECTORS 
OF FLEXIBLE ROTOR-BEARING SYSTEMS WITH SIGLE DISC 


B.l EIGEN VALUES: 

Pi = o 

p i 2i _ l2EI +2n h»h +»h»h+f^h»h»h( mi +»h + m 3 ) + ("V »3 +2m ] n h + m 2 m 3 f 

n h^ 2m x m 1 m l — 


B.2 EIGEN VECTORS: 


U 


with 
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^ U \2 ;/ 22 

1 1 1 


21 


U 22 ~ 
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APPENDIX C 


STATISTICAL ESTIMATION ERRORS 


The statistical errors have been obtained by processing ensembles of data through the parameter estimation 
algorithms. The estimated values of the parameters for a typical ensembles of five signals each have been 
presented here in tabular form. Table C.l concerns the rigid rotor configuration of Chapter 3; Table C.2, C.3 
and C.4 refer to the rotor configurations of Chapters 4, 5 and 6 respectively. 

A typical ensembles of displacement and velocity signals have been shown in Figs C.l and C.2. The 
estimated parameters presented in Table C.l refer to these figures. 

Table C.l Statistical Errors in Single Disc Rigid Rotor Case (Chapter 3) 


Parameters 


co, 


(rads/sec) 2 


Variation 


X (mm' 2 ) 


Variation 


Vertical 5.42 x 10 7 

5.32 x 10 7 
5.37 x 10 7 
5.36 x 10 7 

5,37 x 10 7 

Horizontal 3.21 x 10' 

3.14 x 10 7 
3.11 x 10 7 
3.18 x 10 7 
3.16 x 10 7 


1.85% 


3.11% 


-1.27 x 10 6 
-1.23 x 10 s 
-1.24 x 10 6 
-1.25 x 10 6 
-1.25 x 10 6 
-1.29 x 10 5 
-1.24 x 10 6 
-1.25 x 10 6 
-1.24 x 10 6 
-1.27 x 10 6 


3.15% 


3.88% 


























122 


Table C.2 Statistical Errors in Flexible Single Disc Rotor Case (Chapter 4) 


Parameters 


co 2 (rads/sec)" 


Aj (mm)"' 


a) 2 h (rads/ sec)" 


A, (mm)"' 


Vertical 

Variation 

— 

Horizontal 

Variation 


2.13 x 10 
2.21 x 10 7 

2.14 x 10 : 
2.17 x 10 7 
2.19 x 10 7 


-1.25 x 10° 
-1.24 x 10 6 
-1.21 x 10 6 
-1.25 x 10 s 
-1,23 x 10 6 
1.98 x 10 7 
1.92 x 10 7 

1.96 x 10 7 
1.94 x 10 7 

1.97 x 10 7 
-1.15 x 10 6 
-1.16 x 10 6 
-1.13 x 10 6 
-1.17 x 10 6 
-1,15 x 10 ft 

2.10 

2.06 

2.12 

2.09 

2.11 


3.62% 


3.20% 


3.03% 


3.42% 


2.83% 


2.12 x 10 

2.15 x 10' 
2.17 x 10' 

2.16 x 10' 
2,15 x IQ 7 


-1.41 x 10“ 
-1.44 x 10 s 
-1.42 x 10 6 
-1.43 x 10 6 
-1,42 x 10 6 
2.22 x 10* 

2.24 x 10' 

2.25 x 10 7 
2.21 x 10' 
2.24 x 10' 
-1.36 x 10 6 
-1.32 x 10 6 
-1.35 x 10 6 
-1.34 x 10 6 
-1,32 x 10 6 

2.34 
2.30 

2.35 
2.34 
2.33 


2.30% 


2.08% 


2.94% 


2.13% 
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Table C.3 


Parameters 


k u (N/mm) 


k m (N/mm 3 ) 


lm ,/nj (Kg) 


Statistical Errors in Flexible Multi Disc Rotor Case (Chapter 5) 


Vertical 


1.03 x 10 4 
1.06 x 10 4 

1.04 x 10 4 

1.05 x 10 4 
-3.62 x 10 10 
-3.58 x 10 10 
-3.61 x 10 10 
3.62 x 10 10 
3.59 x IQ 10 

0.21 

0.20 

0.21 

0.21 

0.20 


Variation 


4.76% 


Horizontal 


k L • (N/mm) 

1.04 x 10 4 

1.01 x 10 4 

1.05 x 10 4 

1.03 x 10 4 

1.05 x 10 4 

3.81% 

0.87 x 10 
0.86 x 10 
0.84 x 10 
0.85 x 10 
0.86 x 10 

k NLi (N/mm 3 ) 

-5.10 x 10 10 
-5.08 x 10 10 
-5.06 x 10 10 
-5.09 x 10 10 
-5.06 x 10 10 

0.78% 

-3.50 x 10 
-3.57 x 10 
-3.55 x 10 
-3.57 x 10 
-3.52 x 10 


0.87 x 10 4 
0.84 x 10 4 
0.85 x 10 4 
0.87 x 10 4 
-2.19 x 10 10 
-2.21 x 10 10 
-2.22 x 10 10 
-2.19 x 10 10 
-2,20 x IQ 10 
0.20 
0.21 
0.20 
0.21 
0.20 


3.45% 
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Table C.4 


Statistical Errors in Unbalanced Rotor Case (Chapter 6) 


Speed 

Unbalance 

Parameters estimated "" " ”f 

G) 

m e 

9 

0 


5 

k(x) 

rpm 

gm-cm 

degrees 

degrees 



N/mm 


(at 9 - 0°) 






1800 


2.9 

0.3 

14.23 

0.023 

1.31 xlO 4 -0.92 xlO* u x ; 



2.8 

0.3 

14.42 

0.022 

1.32 xlO 4 -0.91 xlO 10 x 2 


10.5 

2.9 

0.3 

14.32 

0.022 

1.34 xlO 4 - 0.93 x lO'V 

Variation 


2.8 

0.3 

14.30 

0.023 

1.33 x Iff 1 -0.91 x 10 lo xT 



2.9 

0.3 

14.36 

0.022 

1.32 xlO 4 - 0.92 x 10'° x : 



(3.45%) 

(0.00%) 

(1.32%) 

(4.35%) 

(2.24%) (2.15%) 

- 


1.6 

0.1 

20.32 

0.036 

1.41 xlO 4 - 1.10 xlO'V j. 


17.5 

1.6 

0.1 

20.12 

0.036 

1.44 x 10 4 - 1.09 x lO^x 1 



1.6 

0.1 

20.41 

0.035 

1.45 xlO 4 - 1.07 xlO‘°x- 

Variation 


1.6 

0.1 

20.24 

0.036 

1.43 xlO 4 -1.08 xlO 1 V 



1.6 

0.1 

20.31 

0.035 

1.45 xlO 4 - 1.09 xlO'V 



(0.00%) 

(0.00%) 

(1.42%) 

(2.78%) 

(2.76%) (2.73%) 

- 


3.2 

0.2 

27.32 

0.032 

1.60xl0 4 - 1.28xl0‘V 



3.3 

0.2 

27.41 

0.031 

1.62 xlO 4 - 1.27 xlO'V 


24.5 

3.2 

0.2 

27.36 

0.031 

1.63 xlO 4 - 1.25 x 10‘°x- 

Variation 


3.2 

0.2 

27.38 

0.032 

1.61 x 10 4 - 1.26 x 10 10 x : 



3.2 

0.2 

27.40 

0.031 

1.63 x 10 4 - 1.27 x 10 10 x 2 



(3.03%) 

(0.00%) 

(0.33%) 

(3.13%) 

(1.84%) (2.34%) 

1400 


1.7 

0.2 

15.43 

0.063 

0.94x 10 4 -0.94x 10 lo x : 



1.7 

0.2 

15.23 

0.061 

0.91 x 10 4 -0.92x lO^x 1 


10.5 

1.7 

0.2 

15.34 

0.064 

0.93 xlO 4 -0.93 xlO'V 

Variation 


1.7 

0.2 

15.30 

0.062 

0.94X10 4 -0.93 xlO'V 



1.7 

0.2 

15.21 

0.063 

0.93 xlO 4 - 0.92 xlO'V 



(0.00%) 

(0.00%) 

(1.30%) 

(4.69%) 

(3.19%) (2.13%) 

- 


3.2 

0.4 

21.40 

0.063 

0.93 xlO 4 - 1.01 xlO l4 x ; 



3.1 

0.4 

21.32 

0.062 

0.92 xlO 4 - 1.02 xlO'V 


17.5 

3.2 

0.4 

21.37 

0.063 

0.94 x 10 4 - 1.04 x 10 10 x : 

Variation 


3.1 

0.4 

21.41 

0.063 

0.93 x 10 4 - 1.03 x 10'°x : 



3.2 

0.4 

21.35 

0.062 

0.94 xlO 4 - 1.01 xlO'V 



(3.13%) 

(0.00%) 

(0.42%) 

(1.59%) 

(2.13%) (2.88%) 

„ 


1.8 

0.1 

23.56 

0.037 

1.37x 10 J -0.72 x 10 lo x" 



1.8 

0.1 

23.42 

0.036 

1.36 xlO 4 - 0.71 xl0 10 x : 


24.5 

1.8 

0.1 

23.61 

0.037 

1.35 x 10 4 -0.73 x 10 l0 >r 

Variation 


1.8 

0.1 

23.59 

0.036 

1.36 x 10 4 -0.72 x lO 10 ^ 



1.8 

0.1 

23.53 

0.037 

1.37 x 10 4 -0.71 x 10 10 x 2 



(0.00%) 

(0.00%) 

(0.80%) 

j (2.70%) 

(1.46%) (2.74%) 
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